Frage:
Wie simuliere ich NGS-Lesevorgänge und steuere die Sequenzabdeckung?
SmallChess
2017-05-17 21:26:17 UTC
view on stackexchange narkive permalink

Ich habe eine FASTA-Datei mit mehr als 100 Sequenzen wie folgt:

  >Sequence1GTGCCTATTGCTACTAAAA ... >Sequence2GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA ...... Eine Textdatei wie diese:  
  Sequence1 40Sequence2 30 ......  

Ich möchte Paired-End-Lesevorgänge der nächsten Generation für alle simulieren die Sequenzen in meiner FASTA-Datei. Für Sequence1 möchte ich eine 40-fache Abdeckung simulieren. Für Sequence2 möchte ich eine 30-fache Abdeckung simulieren. Mit anderen Worten, ich möchte meine Sequenzabdeckung für jede Sequenz in meiner Simulation steuern.

F: Was ist der einfachste Weg, dies zu tun? Welche Software sollte ich verwenden? Bioconductor?

[Diese Bewertung] (http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html) könnte ein guter Ausgangspunkt sein. Es vergleicht 23 NGS-Simulationswerkzeuge und bietet einen Entscheidungsbaum für die Auswahl eines geeigneten Werkzeugs, das Ihren Anforderungen entspricht.
Welche Leselänge verwenden Sie? Wie lang sind die Sequenzen? Müssen Sie das Abdeckungsziel genau oder mit einiger Wahrscheinlichkeit erreichen?
Klingt nach etwas, das Sie mit dem Biopython-Modul einfach in zehn oder zwanzig Zeilen Python codieren würden.
Ich würde noch ein paar Fragen zu Gregs hinzufügen. Verstehe ich richtig, dass Sie die Sequenzierung der Vorlagensequenz aus der Datei simulieren möchten? Stellen die Sequenzen also Genome dar? Möchten Sie eine Verstärkungsverzerrung in Betracht ziehen? Sequenzierungsfehler? Welche Sequenzierungsplattform möchten Sie simulieren?
Neun antworten:
#1
+6
H. Gourlé
2017-05-30 12:06:03 UTC
view on stackexchange narkive permalink

Ich arbeite an einem Illumina-Sequenzierungssimulator für Metagenomik: InSilicoSeq

Er befindet sich noch in Alpha-Version und ist sehr experimentell, verfügt jedoch über eine Multi-Fasta- und eine Abundance-Datei. Es werden Lesevorgänge aus Ihren Eingabegenomen mit unterschiedlichen Abdeckungen generiert.

Aus der Dokumentation:

  iss generate --genomes genomes.fasta --abundance reichlich_file.txt \ - model_file HiSeq2500 - Ausgabe von HiSeq_reads  

Wobei:

  # multi-fasta file>genome_AATGC ... >genome_BCCGT ...... # Abundanzdatei (Gesamthäufigkeit) muss 1 sein!) Genom_A 0.2genom_B 0.4 ...  

Ich habe es nicht für die Abdeckung entworfen, sondern für die Fülle des Genoms in einem Metagenom, daher müssen Sie möglicherweise a kleines bisschen Mathe;)

#2
+4
Ian Sudbery
2017-05-18 13:57:06 UTC
view on stackexchange narkive permalink

Das Polyester -Bio-Leiterpaket kann dies. Es heißt, es simuliert RNA-seq-Lesevorgänge, aber ich weiß nicht, ob sich das wirklich von anderen NGS-Lesevorgängen unterscheidet.

Es kann eine Reihe von Fehler- und Bias-Modellen verwenden oder sie aus einem Datensatz lernen.

Könnten Sie das etwas näher erläutern? Wie würde das OP es für seine Dateien verwenden? Welcher Bereich von Fehlern und Verzerrungen? Wie würden wir sie verwenden? Was wäre eine vernünftige Standardeinstellung? Vielleicht wäre ein minimales Arbeitsbeispiel hilfreich.
#3
+4
Kamil S Jaron
2017-05-18 22:32:11 UTC
view on stackexchange narkive permalink

Dieses Python-Skript verwendet eine Fasta-Datei und eine Tsv-Datei mit Zählwerten und druckt die Sequenzen in den Fasta-Dateien so oft aus, wie es in der Tsv-Datei angegeben ist (unter der Annahme des in der Frage angegebenen Formats). Wenn also bar.tsv und foo.fasta ​​code> Ihre Dateien sind:

  aus Bio-Import SeqIOrepeat = {} für Zeile in open ( "bar.tsv"): seq_id, Coverage = line.split () repeat [seq_id] = int (Coverage) für seq_record in SeqIO.parse (foo.fasta, "fasta"): für i in range (repeat.get ( seq_record.name, 0)): print (">", seq_record.name, "_", i, sep = '') print (seq_record.seq)  
Können Sie erklären, was dies bewirkt und wie es funktioniert? Es sieht so aus, als würden Sie nur die Originalsequenz drucken, und das kann nicht richtig sein. Das OP muss Lesevorgänge simulieren und nicht nur die gleiche Sequenz N-mal wiederholen.
Autsch, dann habe ich die Frage falsch verstanden. Ich denke, dass die Sequenzen die Lesevorgänge sind, die simuliert werden sollten und 40x seitlich ist. Ich möchte diese Sequenz 40 Mal haben. Ich habe die Antwort bearbeitet, um zumindest Verwirrung zu vermeiden.
Ich weiß nicht, ob du es falsch verstanden hast oder ich. Aber einer von uns hat es getan :) Danke für die Bearbeitung, zumindest jetzt kann das OP sehen, was dies tut und selbst entscheiden.
#4
+4
Daniel Standage
2017-05-23 00:13:50 UTC
view on stackexchange narkive permalink

Das wgsim-Paket von Heng Li (von BWA und samtools) ist mein Werkzeug zur Simulation von Illumina-Lesevorgängen. Es bietet keine bequeme Möglichkeit, die unterschiedliche Abdeckung über verschiedene Sequenzen hinweg zu simulieren, aber es sollte nicht zu schwierig sein, wgsim mehrmals auszuführen und den gewünschten Abdeckungsgrad für jede interessierende Sequenz zu generieren.

I. würde ein Python-Skript implementieren, um Ihre Testdatei zu schlürfen, und wgsim (unter Verwendung des Moduls subprocess ) für jede Sequenz aufrufen. Dies erfordert wahrscheinlich, dass Sie jede Sequenz in einer separaten Datei haben. :-(

Könnten Sie dies erweitern, indem Sie einen Beispielbefehl hinzufügen, der zeigt, wie eine Fasta-Eingabedatei verwendet und das Fastq erstellt wird? Dies fühlt sich eher nach (hilfreichen) Anweisungen an, wo eine Antwort gefunden werden kann, als nach einer Antwort.
#5
+2
Gabriel Renaud
2017-05-17 21:36:37 UTC
view on stackexchange narkive permalink

Mir ist keine Software bekannt, die dies direkt ausführen kann, aber ich würde die Fasta-Datei in eine Sequenz pro Datei aufteilen, sie in BASH durchlaufen und ART den Sequenzsimulator (oder eine andere) für jede Sequenz aufrufen.

Könnten Sie das näher erläutern? Dies ist in diesem Zustand wirklich keine sehr nützliche Antwort. Was ist Kunst? Wo kann ich es finden? Wie benutze ich es? Wie funktioniert iot? Führt es Fehler ein? Können wir die Fehlerrate kontrollieren? Nimmt es wirklich keine Multifasta-Dateien? Wenn nicht, wird das Schleifen über mehrere tausend Dateien in Bash sehr, sehr langsam sein. Es kann bessere Wege geben.
> Was ist KUNST? Es ist schwer zu definieren, ich würde sagen, es handelt sich um verschiedene Aktivitäten im Zusammenhang mit der Herstellung von Werken visueller oder auditorischer Produkte, die darauf abzielen, das Innere des Autors zu externalisieren ... oh, Sie meinen den Sequenzsimulator? Das Papier finden Sie hier: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/> Wo finde ich es? siehe Papier> Wie verwende ich es? Handbuch lesen> Wie funktioniert iot? Papier lesen.> Führt es Fehler ein? Papier / Handbuch lesen Können wir die Fehlerrate kontrollieren? Papier lesen.
> Nimmt es wirklich keine Multifasta-Dateien? Wenn nicht, wird das Schleifen über mehrere tausend Dateien in Bash sehr, sehr langsam sein. Es kann bessere Möglichkeiten gebenWo hat jemand angegeben, dass Multifasta nicht benötigt wird? OP wollte eine sehr spezifische Abdeckung pro Fasta-Datensatz.
Nein, ich wollte es mir nicht erklären. Ich meinte, bitte bearbeiten Sie Ihre Antwort, damit sie diese Informationen enthält. Wir versuchen, Antworten so umfassend wie möglich und schriftlich zu haben, damit sie ohne externe Referenzen auf den Stack Exchange-Sites allein stehen können. Kommentare sind schwer zu lesen, leicht zu übersehen und führen zu einem schrecklichen Durcheinander von Biostars, bei denen Sie keine Ahnung haben, wo die eigentliche Antwort liegt. Bitte bearbeiten Sie Ihre Frage zur Klärung.
Deine Antwort, meinte ich. Es tut uns leid.
#6
+2
Greg
2017-05-17 22:55:25 UTC
view on stackexchange narkive permalink

Wenn Sie mit einer gewissen Zufälligkeit einverstanden sind, können Sie mit einer Poisson-Zufallsvariablen Lesevorgänge aus Ihrer Sequenzdatei generieren. Sie müssen einige Berechnungen durchführen, um herauszufinden, welcher Lambda-Wert verwendet werden soll, damit die erwartete Abdeckung bei jedem Basenpaar in Ihrem Lesevorgang mit der in Ihrer Textdatei festgelegten übereinstimmt.

Zum Beispiel Sie haben eine Sequenz S mit einer Länge von 1.000, eine Leselänge von 50 und eine Insertgröße von 100. Für jede Basis b in S wird eine Poisson-Zufallsvariable p erzeugt. Sie generieren dann p Lesevorgänge von Basis b bis b + 50. Generieren Sie dann den gepaarten Lesevorgang ab b + 50 + 100.

Auch hier müssten Sie damit spielen, um herauszufinden, welches Lambda verwendet werden soll. Dies würde Ihnen jedoch im Grunde das geben, was Sie möchten, solange Sie nicht genau die Abdeckung haben, auf die Sie abzielen für jeden Lesevorgang.

Und es wird auch realistischer sein, den Einfügungsabstand von der Poisson-ähnlichen Verteilung zu simulieren! Die Sache ist, dass dies eine Neuerfindung eines Fahrrads ist. Es gibt viele Simulationswerkzeuge.
Wäre das Ziehen aus einer Normalverteilung für die Einfügungsgröße nicht besser? Und ja, ich weiß, dass es viele Tools gibt, aber dies wäre eine Möglichkeit, es selbst zu tun, wenn Sie interessiert wären.
Nach meiner Erfahrung ist normal etwas weniger flexibel für die Anpassung - die Dichte ist verzerrt, z. B. http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png. Wahrscheinlich ist ein negatives Binom eine gute Wahl, aber ich habe nie versucht, es an die tatsächliche Dichte der Insertgröße anzupassen ...
#7
+1
Karel Brinda
2017-06-21 19:46:49 UTC
view on stackexchange narkive permalink

Mit RNFtools (ab Version 0.3.1) ist es jetzt einfach, NGS-Lesevorgänge zu simulieren und gleichzeitig die Sequenzabdeckung zu steuern. Weitere Informationen finden Sie im -Tutorial, insbesondere im Abschnitt Sequenzextraktion.

Vorbereitung der Umgebung

Installieren Sie zunächst BioConda und fügen Sie die erforderlichen Kanäle hinzu. Installieren Sie dann entweder RNFtools in der Standard-Conda-Umgebung

  conda install rnftools  

oder erstellen und aktivieren Sie eine separate Conda Umgebung (vorzuziehen)

  conda create -n rnftools rnftoolssource aktiviere rnftools  

Simulation

Angenommen, Sie haben eine Referenzdatei ref.fa und eine durch Tabulatoren getrennte Coverage-Datei Coverage.tsv (z. B. die aus Ihrem Beispiel). . Dann erledigen die folgenden RNFtools Snakefile den gewünschten Job:

  importieren Sie rnftoolsimport csvrnftools.mishmash.sample ("simulations_mit_abdeckungssteuerung") , read_in_tuple = 1) fa = "ref.fa" tsv = "Coverage.tsv" mit open (tsv) als f: table = csv.reader (f, delimiter = '\ t') für seqname, cov in table: rnftools .mishmash.DwgSim (fasta = fa, Sequenzen = [seqname], Coverage = float (cov), read_length_1 = 10, # Schnelltest mit supershort read read_length_2 = 0,) include: rnftools.include () Regel: Eingabe: rnftools. input ()  

Wenn Sie diese Datei ( Snakefile ) speichern und snakemake ausführen, simuliert RNFtools Lesevorgänge mit DWGsim mit den definierten Abdeckungen Ihre Textdatei und speichern Sie alle simulierten Lesevorgänge in simulation_with_coverage_control.fq .

Sie können mit allen Parametern spielen. Insbesondere können Sie einen anderen Simulator verwenden (z. B. Art-Illumina mit rnftools.mishmash.ArtIllumina ). Weitere Informationen finden Sie in der RNFtools-Dokumentation.

#8
  0
winni2k
2017-06-12 01:27:23 UTC
view on stackexchange narkive permalink

Ein weiteres Lesesimulationstool der nächsten Generation ist gemsim. Ich habe es nicht getestet, aber es würde mich interessieren, ob jemand Erfahrung damit hat.

#9
-1
Karel Brinda
2017-05-30 00:46:25 UTC
view on stackexchange narkive permalink

Sie können Ihre FASTA-Datei nacheinander aufteilen, indem Sie

  split -a 6 -p '^ >' your_file.fa seq_  

verwenden und verwenden Sie dann einen vorhandenen Lesesimulator, der die Abdeckung unterstützt (ART, DWGsim usw.). Wenn Sie möchten, dass alle Lesevorgänge gemischt werden (nicht nach der ursprünglichen Reihenfolge sortiert), können Sie RNFtools verwenden.

Bearbeiten 1:

Als @terdon Der vorherige Befehl funktioniert nur unter OS X. Ein analoger Einzeiler für Linux (jedoch mit einem etwas anderen Namensschema, bei dem Zahlen anstelle von Buchstaben verwendet werden) kann

  csplit -f seq_ -n 6 your_file.fa '/ ^ > /' {* sein }  

Damit dieser Befehl auch unter OS X funktioniert, müssen Coreutils installiert werden (z. B. mit Brew) und dann gcsplit anstelle von csplit verwendet werden .

Bearbeiten 2:

Sobald FASTA durch Sequenzen aufgeteilt ist, wird die Simulation unkompliziert und es können viele verschiedene Ansätze verwendet werden. Mein Favorit ist GNU Parallel. Stellen Sie sich vor, Sie haben Ihre Coverages in einer Textdatei mit dem Namen covs.txt in separaten Zeilen und in derselben Reihenfolge wie die Sequenzen in your_file.fa , z. B.

  4030...  

Anschließend können Sie mit DWGsim Lesevorgänge aus den Originalsequenzen simulieren, indem Sie

  ls -1 seq_ * | füge covs.txt ein - \ | parallel -v --colsep '\ t' dwgsim -1 100 -2 100 -C {1} {2} sim_ {2}  

und füge die erhaltenen FASTQ-Dateien zusammen mit:

  cat sim_seq _ *. bwa.read1.fastq > reads.1.fqcat sim_seq _ *. bwa.read2.fastq > reads.2.fq  

Eine Möglichkeit Die Gefahr dieses Ansatzes besteht darin, dass wir angenommen haben, dass die Anzahl der seq_ * -Dateien der Anzahl der Zeilen in covs.txt entspricht, was möglicherweise nicht der Fall ist (aus Versehen). Wir sollten dies vor dem Simulationsschritt überprüfen, z. B.:

  [["$ (ls -1 seq_ * | wc -l)" == "$ (cat covs.txt | wc -l) "]] \ || echo "Inkorrekte Anzahl von Zeilen in covs.txt"  

Eine weitere Einschränkung ist, dass die simulierten Lesevorgänge nicht in zufälliger Reihenfolge erfolgen (sie sind nach Quellsequenzen gruppiert).

Ich fürchte, Sie beantworten die Frage nicht. Zunächst verwenden Sie nicht standardmäßige Sortieroptionen. Ich nehme an, Sie verwenden einen Mac? Das `-p` ist eine reine BSD-Option und auf anderen * nix-Systemen AFAIK nicht verfügbar. Das heißt, während dies wahrscheinlich funktioniert, um die Datei auf einem MacOS-System zu teilen, ist das Teilen trivial und es gibt viele, viele Möglichkeiten, dies zu tun. Das Schwierige ist, die Lesevorgänge zu simulieren, und Sie erklären nicht wirklich, wie das OP das tun kann.
Danke für die Bearbeitung. Das Hauptproblem hierbei ist jedoch, dass Ihre Antwort die gestellte Frage nicht beantwortet. Sie antworten: "Wie kann ich eine Multifasta-Datei in viele Einzelsequenzdateien aufteilen?", Aber die Frage betrifft die Simulation von Lesevorgängen. Das Teilen der Datei kann * Teil * einer Antwort sein oder auch nicht, aber es ist sicherlich keine Antwort.
Oh, und für einen tragbaren Weg, ohne eines der vielen Spezialwerkzeuge zu verwenden, können Sie "awk" verwenden. So etwas wie "awk -vRS ="> "-F" \ n "NR> 1 {print"> "$ 0 >> $ 1" .fa "}" file.fa ". Ihr "csplit" (und Ihr "split") Ansatz ist jedoch in der Tat einfacher.


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