Frage:
Wie kann eine Teilmenge von bam sicher und effizient in fastq konvertiert werden?
Kamil S Jaron
2017-07-26 16:54:47 UTC
view on stackexchange narkive permalink

Frage

Wie kann ich Lesevorgänge aus einer bam -Datei (erstellt von bwa-mem ) nach fastq extrahieren? eine Liste von Referenzsequenzen zum Herausfiltern gegeben?

Mögliche Schwierigkeiten

  • Beibehalten der FR-Ausrichtung von Lesevorgängen am Paarende (in bam sind alle Sequenzen Referenzsequenzen)
  • Halten von R1 und R2 liest
  • Halten von Qualitätsbewertungen in derselben Codierung wie Original-Fastq (in meinem Fall Standard-Illumina-Phrasen-Bewertungen)
  • bam can be (ana ist normalerweise) nach Koordinaten sortiert

Fast gibt es Lösungen

  1. Sujais Perl-Lösung in Die Blobologie macht genau das Gegenteil - sie erhält Lesungen aus der Referenzliste (also könnte ich die Liste einfach umkehren). Der Nachteil ist, dass das Skript eine verschachtelte fq -Datei ausgibt. erfordert eindeutige Namen von Partnern, andernfalls gehen R1 / R2-Informationen verloren.

  2. samtools + grep sie alle aus fastq-Dateien

  3. ol>

    Erstellen Sie eine Liste der gelesenen Namen, die nicht gefilterten Gerüsten zugeordnet sind. (cut extrahiert nur gelesene Namen, uniq reduziert die gelesenen Namen von Paarende, wenn sie gleich sind). Dann grep grep liest Namen aus Fastq-Dateien und entfernt - Trennzeichen zwischen Treffern

      samtools view foo.bam | grep -vf list_of_scaffols_filter \ | cut -f 1 | uniq > list_of_reads_to_keepgrep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^ - $" > foo_R1_filtered_bash.fqgrep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^ - $" > foo_R1_filtered_bash.fq  
    1. filter bam & picard tools
    2. ol>

      Oder ich könnte nur den Filterteil machen und Picard-Tools (Picard.SamToFastq) verwenden, aber wie üblich vermeide ich Java so weit wie möglich. Ich denke

        samtools view foo.bam | grep -vf list_of_scaffols_filter \ | java -jar picard.jar SamToFastq INPUT = / dev / stdin \ FASTQ = foo_R1_filtered_bash.fq SECOND_END_FASTQ = foo_R2_filtered_bash.fq  

      Die erste Lösung funktioniert bei mir nicht wirklich, da ich nicht alle Lesevorgänge in der BAM-Datei umbenennen und die R1 / R2-Informationen beibehalten möchte (da R1 und R2 unterschiedliche Fehlerprofile haben). Beide Lösungen 2 und 3 finde ich etwas ungeschickt und ich bin mir nicht sicher, ob sie allgemein sind. Ich kann ein unerwartetes Verhalten bekommen, wenn einer der Lesevorgänge nicht als zweiter zugeordnet wird. Beide leiten denselben Filterungsschritt weiter.

      Ich habe mich über eine Pysamlösung gewundert. Ich denke, es wird viel langsamer sein, aber zumindest wird es viel klarer und vielleicht allgemeiner. So etwas wie in Bam-Datei in Fasta-Datei konvertieren - es gibt eine Pysam-Lösung für Fasta (nicht Fastq), fast da ...

      Story

      Ich habe sehr fragmentiertes Referenzgenom. Einige der Gerüste sind zu klein, um damit zu arbeiten, und einige von ihnen sind Verunreinigungen (identifiziert mit blobtools). Ich möchte Lesevorgänge trennen, die verschiedenen Gruppen zugeordnet sind, um Verunreinigungen, kurze Gerüste und Gerüste zu trennen, die für die nachgeschaltete Analyse verwendet werden. Der Grund dafür ist, dass, wenn wir alle Lesevorgänge einer gefilterten Referenz (0,7 - 0,8 des ursprünglichen Genoms) zuordnen, die meisten von ihnen (0,95 - 0,99) immer noch einen Ort finden, an dem sie abgebildet werden. Daher gibt es 0,2 - 0,3 falsch platzierte Lesevorgänge wird offensichtlich die Downstream-Analyse beeinflussen müssen, wie das Aufrufen von Varianten.

      Diese Filteridee basiert auf der Logik, dass wenn die gefilterte duplizierte Genomregion einige kleine Unterschiede enthält, sie ihre Lesevorgänge anziehen (und wenn ich sie filtere, ich verbessert das Aufrufen von Varianten) und wenn sie genau gleich sind, werden ihnen zufällig Lesevorgänge zugewiesen, sodass dies keinen Schaden anrichtet.

Was möchten Sie tun, wenn nur ein Lesevorgang einem Paar einem Gerüst / Contig zugeordnet ist, das Sie ausschließen möchten? Ehrlich gesagt würde ich einfach die BAM-Datei sortieren und ein bisschen Python schreiben, um die Konvertierung und Filterung durchzuführen, aber ich bin nur faul.
Ich bin mir noch nicht sicher, aber ich denke, es ist besser, alle Paare zu behalten, bei denen mindestens eines der Paare einem nicht gefilterten Gerüst zugeordnet ist. - Wenn R1 / 2 dieselben Namen hat, ist Lösung 2 genau so werde mich benehmen. Wenn die Namen unterschiedlich sind, filtere ich die unterschiedliche Anzahl von R1- und R2-Lesevorgängen.
Zwei antworten:
#1
+5
Devon Ryan
2017-07-26 18:13:32 UTC
view on stackexchange narkive permalink

Mir ist kein vorgefertigtes Programm dafür bekannt, deshalb habe ich eines für Sie geschrieben. Dies nimmt eine BAM-Datei mit jeder Bestellung und erzeugt ordnungsgemäß geordnete komprimierte Fastq-Dateien mit der von Ihnen angeforderten Filterung. Intern durchläuft dies alle Einträge in der BAM-Datei (ignoriert sekundäre / ergänzende Einträge und diejenigen, bei denen beide Partner Ihrer Filterliste zugeordnet sind), speichert die ordnungsgemäß ausgerichtete Sequenz / Qualität / den gelesenen Namen in einem Puffer und speichert diesen Puffer dann Eintrag auf die Festplatte, sobald der Partner gefunden wurde. Dies sollte einigermaßen performant sein (hey, es ist Python, um nicht zu viel zu erwarten). Wenn Sie jedoch BAM-Dateien indiziert haben, können Sie sich überlegen, wie Sie dies schneller ausführen können.

Überprüfen Sie dies die Ausgabe, da ich nur einen Test ausgeführt habe.

Sie waren 3 Minuten schneller und Ihr Code sieht viel besser aus. Ich werde es ausprobieren ... Vielen Dank.
Der einzige wirkliche Vorteil meines Codes besteht darin, dass er nicht in der BAM-Datei herumspringt, um nach Partnern zu suchen, und dass die BAM-Datei nicht sortiert werden muss. Wenn Sie Dateien sortiert haben, ist es meiner Meinung nach schneller, eine Liste von Contigs zu erstellen, die nicht ausgeschlossen sind, fügen Sie "*" hinzu und iterieren Sie dann wie in meinem Programm darüber. Das würde mehr Zeit sparen.
#2
+2
Kamil S Jaron
2017-07-26 18:17:23 UTC
view on stackexchange narkive permalink

Ok, ich habe einen etwas brutalen Pysam / BioPython-Parser geschrieben, der den Index von bam verwendet, um die richtige Reihenfolge der Lesepaare für R1 / R2-Dateien und das bitweise Flag zu erhalten. Es sollte nicht allzu schwierig sein, jetzt komplexere Filterregeln hinzuzufügen.

  #! / Usr / bin / env python3 # 1. arg - indizierte BAM-Datei # 2. arg - Liste der Header Filter # 3. arg - Namensmuster für die Ausgabe readimport osimport sysimport pysamfrom Bioimport SeqIO, Seq, SeqRecordsamfile = pysam.AlignmentFile (sys.argv [1], "rb") header_set = set (line.strip () für line in open (sys.argv [2])) base = sys.argv [3] out_R1 = base + 'R1_filtered.fq'out_R2 = base +' R2_filtered.fq'with open (out_R1, mode = 'w') als R1, öffne (out_R2, mode = 'w') als R2: für den Eintrag in der Samfile: wenn entry.is_read1 und nicht entry.reference_name in header_set: # pait pair entry_R2 = samfile.mate (entry) # mache eine Sequenzgymnastik mit R1 seq_R1 = Seq.Seq (entry.seq) wenn entry.is_reverse: seq_R1 = seq_R1.reverse_complement () # mache eine Sequenzgymnastik mit R2 seq_R2 = Seq.Seq (entry_R2.seq) wenn entry_R2.is_reverse : seq_R2 = seq_R2.reverse_complement () R1.write ('@' + entry.qname + '\ n' + str (seq_R1) + '\ n + \ n' + entry.qqual + '\ n') R2.write ( '@' + entry_R2.qname + '\ n' + str (seq_R2) + '\ n + \ n' + entry_R2.qqual + '\ n') samfile.close ()  

Es gibt einige hässliche Teile, bitte zögern Sie nicht, es schöner zu machen.



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