Frage:
Abrufen eindeutig zugeordneter Lesevorgänge aus der BWA-Mem-Ausrichtung
gringer
2017-06-07 03:23:31 UTC
view on stackexchange narkive permalink

Dies basiert auf einer Frage von betsy.s.collins zu BioStars. Den Originalbeitrag finden Sie hier.

Hat jemand Vorschläge für andere Tags oder Filterschritte für BWA-generierte BAM-Dateien, die nur zum Lesen verwendet werden können? Karte zu einem Ort? Eine Beispielanwendung wäre das Finden von Seeds für den TULIP Assembler / Scaffolder, der am besten für Lesevorgänge geeignet ist, die eindeutigen genomischen Positionen zugeordnet sind.

Das Konzept der "eindeutig zugeordneten Lesevorgänge" lautet ein geladener Begriff, und die meisten Quellen schlagen vor, dass das Filtern nach MAPQ den Trick tun sollte. Dieser Ansatz scheint jedoch nicht zu funktionieren, wenn BWA als Lese-Mapper verwendet wird.

Ich bin ein bisschen frech mit diesem Beitrag, hauptsächlich um zu beurteilen, ob die Community dies für eine gute Idee hält oder nicht. Dies ist eine Frage, die gerade von einem völlig neuen Benutzer auf BioStars gestellt wurde. Die Frage ist sehr spezifisch und enthält eine Menge guter Geschichten. Ich denke, sie passt gut zum Stack Exchange-Modell.
Eindeutig zugeordnete Lesevorgänge (d. H. Ein Lesevorgang, der einem einzelnen Ort im Genom zugeordnet ist) werden manchmal bevorzugt, wenn Analysen ausgeführt werden, die von der Quantifizierung der Lesevorgänge abhängen und nicht nur von der Abdeckung (z. B. RNASeq). Doppelte Lesevorgänge erfordern zusätzliche Filterung oder Normalisierung auf der Analyseseite, und die meisten nachgeschalteten Programme berücksichtigen keine Konzepte wie "ein Fünftel der gelesenen Karten irgendwo * hier * auf Chromosom 5 und vier Fünftel der gelesenen Karten irgendwo * hier * auf * Chromosom 14. "
Der Fragesteller bezieht sich auf genomisch eindeutige Alignments, die sich von mehreren möglichen Alignments an einem einzelnen Ort unterscheiden. Wenn sich das einzige Mehrfach-Alignment am selben genomischen Ort befindet, sind diese Alignments im genomischen Sinne immer noch eindeutig. Die Unmöglichkeit, die genaue, korrekte Ausrichtung zu finden, ist ein bekanntes Problem, und es gibt einige nachgeschaltete Methoden (z. B. Linksnormalisierung), um sicherzustellen, dass mehrere lokale Ausrichtungen und / oder Sequenzierungsfehler eine geringere Auswirkung haben.
Vielen Dank. Wenn Sie Text von einer anderen Stelle wörtlich reproduzieren möchten, sollten Sie ihn immer in einen Anführungszeichen setzen (verwenden Sie den Editor oder fügen Sie am Anfang der Zeile ein `>` hinzu), um deutlich anzuzeigen, dass er in Anführungszeichen steht. Wenn nicht, besteht ein hohes Risiko, dass Ihr Beitrag wegen Plagiats gelöscht wird, selbst wenn Sie (wie Sie) die Quelle angeben.
Das Originalplakat hat nun auf Biostars geantwortet und erwähnt, dass sie mit der MAPQ-Filterung beginnen wird (zusammen mit dem Entfernen zusätzlicher Ausrichtungen).
Sie erwähnte auch, dass ihr * tatsächliches * Problem darin bestand, dass Geneious- und BWA-Zuordnungen nicht gleich waren (was verständlich war, wenn jemand in die "Textwand" las, die gepostet wurde); Selbst das Ausschließen von Lesevorgängen mit mehreren Zuordnungen konnte dies nicht beheben.
Zwei antworten:
#1
+9
gringer
2017-06-07 09:21:29 UTC
view on stackexchange narkive permalink

Um alle möglichen mehrfach zugeordneten Lesevorgänge aus einer BWA-zugeordneten BAM-Datei auszuschließen, müssen Sie anscheinend grep für die nicht komprimierten SAM-Felder verwenden:

  samtools view -h mapped.bam | grep -v -e 'XA: Z:' -e 'SA: Z:' | samtools view -b > unique_mapped.bam  

Erklärung folgt ...

Ich gehe von einer Situation aus, in der einem Bioinformatiker eine zugeordnete BAM-Datei präsentiert wird produziert von BWA und hat keine Möglichkeit, die ursprünglichen Lesungen zu erhalten. Eine aufwändige Lösung wäre, die zugeordneten Lesevorgänge aus der BAM-Datei zu extrahieren und mit einem anderen Mapper neu zuzuordnen, der den MAPQ-Score verwendet, um mehrere Zuordnungen anzuzeigen.

... aber was wäre, wenn dies nicht der Fall wäre? möglich?

Mein Verständnis der Ausgabe von BWA ist, dass ein Lesevorgang, der perfekt auf mehrere genomische Orte abgebildet wird, für beide Orte einen MAPQ-Wert (High Mapping Quality) erhält. Viele Leute erwarten, dass ein Lesevorgang, der mindestens zwei Orten zugeordnet ist, (bestenfalls) eine 50% ige Wahrscheinlichkeit hat, einem dieser Orte zugeordnet zu werden (d. H. MAPQ = 3). Da BWA dies nicht tut, ist es schwierig, mehrfach zugeordnete Lesevorgänge aus BWA-Ergebnissen mithilfe des MAPQ-Filters herauszufiltern, der für andere Aligner funktioniert. Dies ist wahrscheinlich der Grund, warum die aktuelle Antwort auf Biostars [ samtools view -bq 1 ] wahrscheinlich nicht funktioniert.

Hier ist eine Beispielzeile aus einer BWA-Mem-Ausrichtung, die ich habe gerade gemacht. Dies sind Illumina zu einem Parasiten Genom kartiert liest, die eine hat viel von sich wiederholenden Sequenz:

  ERR063640.7 16 tig00019544 79.974 21 21M2I56M1I20M * 0 0 TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D NM: i : 4 MD: Z: 83A13 AS: i: 77 XS: i: 67 XA: Z: tig00019544, -78808,21M2I56M1I20M, 6; tig00019544, -84624,79M1I20M, 6; tig00019544, -79312,33M4I42M1I20M, 8; Code> 

BWA mem hat festgestellt, dass dieser bestimmte Lesevorgang ERR063460.7 mindestens drei verschiedenen Speicherorten zugeordnet ist: tig00019544, tig00019544 und tig00019544. Beachten Sie, dass der MAPQ für diesen Lesevorgang 21 ist. Obwohl der Lesevorgang mehreren Speicherorten zugeordnet ist, kann MAPQ nicht verwendet werden, um dies zu bestimmen.

Die alternativen Speicherorte werden jedoch durch das Vorhandensein des XA -Tag im Abschnitt "Benutzerdefinierte Felder" der SAM-Ausgabe. Vielleicht kann nur das Filtern nach Zeilen, die das XA -Tag enthalten, mehrfach zugeordnete Lesevorgänge ausschließen. Die Manpage samtools view schlägt vor, dass -x ein bestimmtes Tag herausfiltert:

  $ samtools view -x XA output.bam | grep '^ ERR063640 \ .7 [[: space:]]' ERR063640.7 tig00019544 16 79.974 21 21M2I56M1I20M * 0 0 TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D NM: i: 4 MD: Z: 83A13 AS: i: 77 XS: i: 67  

... so dass das Tag herausgefiltert wurde (dh das Tag ist in der SAM-Ausgabe nicht mehr vorhanden), aber nicht das read . Das FLAG-Feld enthält keine nützlichen Bits, um mehrere genomische Zuordnungen anzuzeigen (von denen ich weiß, dass sie gefiltert werden können, um auch das Lesen auszuschließen). Daher muss ich auf andere Maßnahmen zurückgreifen.

In diesem speziellen Fall Ich kann grep -v für die unkomprimierte SAM-Ausgabe verwenden, um Ausrichtungslinien mit dem Tag XA auszuschließen (und anschließend erneut auf BAM zu komprimieren, um aufgeräumt zu sein):

  $ samtools view -h output.bam | grep -v 'XA: Z:' | samtools view -b > output_filtered.bam $ samtools view output_filtered.bam | grep '^ ERR063640 \ .7 [[: space:]]' <no output>  

Hurra! liest gefiltert. Abgesehen davon hat diese grep -Suche eine ziemlich erhebliche Rechenlast: Sie sucht irgendwo in der Zeile nach einer Zeichenfolge mit dem Text XA: Z: und tut dies nicht tatsächlich jede Situation erfassen. Einige masochistische Personen könnten zu einem späteren Zeitpunkt vorbeikommen und entscheiden, dass sie alle ihre Lesevorgänge HAXXA: Z: AWESOME !: <readNumber> aufrufen. In diesem Fall wäre eine Optimierung dieser grep-Suche erforderlich um sicherzustellen, dass vor dem XA: Z: ein Leerzeichen (oder genauer gesagt ein Tabulatorzeichen) steht. Jetzt überprüfe ich any doppelte gelesene Namen, nur um sicherzugehen:

  $ samtools view output_filtered.bam | awk '{print $ 1}' | sortieren | uniq -dERR063640.1194ERR063640.1429ERR063640.1761ERR063640.2336ERR063640.2825ERR063640.3458ERR063640.4421ERR063640.4474ERR063640.4888ERR063640.49ERR063640.4974ERR063640.5070ERR063640.5130ERR063640.5300ERR063640.5868ERR063640.6116ERR063640.6198ERR063640.6468ERR063640.6717ERR063640.6797ERR063640.7322ERR063640.750ERR063640.7570ERR063640. 7900ERR063640.8115ERR063640.8405ERR063640.911ERR063640.9206ERR063640.9765ERR063640.9986  

Oh ... verdammt. Ich frage mich, was sie sind:

  $ samtools view output_filtered.bam | grep '^ ERR063640.3458 [[: space:]]' ERR063640.3458 tig00002961 16 5402 60 0 0 * 58S38M AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA DD6HFFJBKFH = KDILKLGLJEKLKGFJIH8IKHLLMJEK: L: HBGJIHJKFLLKIHJDHLNKCK; KMKGMFKJILIIIMKI9JLKKHEJFII CC NM:? i: 0 MD: Z: 38 AS: i : 38 XS: i: 0 SA: Z: tig00002377,202353, -, 14M3I5M1I35M38S, 19,5; ERR063640.3458 2064 tig00002377 202353 14M3I5M1I35M38H 19 * 0 0 AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA DD6HFFJBKFH = KDILKLGLJEKLKGFJIH8IKHLLMJEK: L: HBGJIHJKFLLKIHJ NM: i: 5 MD: Z : 5G48 AS: i: 35 XS: i: 27 SA: Z: tig00002961,5402, -, 58S38M, 60,0;  

Aha! Zusätzliche Ausrichtungen, die das offizielle SA -Tag verwenden [andere kanonische Ausrichtungen in einer chimären Ausrichtung]. Dies scheinen Situationen zu sein, in denen ein einzelner Lesevorgang aufgeteilt und mehreren Speicherorten zugeordnet wurde. Beachten Sie, dass in diesem Fall beide Ausrichtungen immer noch MAPQ-Werte von über 3 aufweisen. Es hört sich so an, als würde der Fragesteller auch diese Situationen loswerden wollen. Dieses Mal gibt es auch Standard-Flag-Felder, um diese Situationen zu behandeln (0x800: sekundäre Ausrichtung). Es reicht jedoch nicht aus, nur die zusätzliche Ausrichtung zu filtern, da beide Lesezuordnungen entfernt werden sollten und nicht nur die (die), die zufällig als sekundär gekennzeichnet wurden.

Glücklicherweise scheint BWA das SA -Tag in alle Lesevorgänge einzufügen, die zusätzliche Ausrichtungen enthalten (wenn dies nicht der Fall ist, bin ich sicher, dass mich jemand korrigieren wird). Daher füge ich in der SA -Suche einen zusätzlichen grep -Filter hinzu:

  $ samtools view -h output.bam | grep -v -e 'XA: Z:' -e 'SA: Z:' | samtools view -b > output_filtered.bam $ samtools view output_filtered.bam | awk '{print $ 1}' | sortieren | uniq -d<no output>  

Fertig. Kinderleicht! </sarcasm>

... diese ursprüngliche "mühsame" Lösung, einen anderen Aligner zu verwenden, sieht jetzt nicht so schlecht aus.

Das ist eine ordentliche Lösung und gut argumentiert. Ich denke nicht, dass es so aufwändig ist, wenn Sie das Ganze auch auf ein paar Zeilen Shell-Skript zusammenfassen. :) :)
"SA" bezeichnet eher zusätzliche Ausrichtungen als sekundäre Ausrichtungen.
Die Verwendung eines MAPQ zum Herausfiltern von Multimappern besteht darin, dass davon ausgegangen wird, dass ein Multimapper "gleich gut zugeordnet" ist, was selbst für die hier vorgestellten Beispiele nicht der Fall ist. Ich habe Ihre Antwort als akzeptierte Antwort auf Biostars festgelegt, aber dort gebeten, ein wenig zu klären, was genau gewünscht wird.
Ja, es ist nicht schwer, wenn Sie wissen, was Sie erwartet. Deshalb habe ich versucht, alle Probleme zu lösen, die auftreten könnten.
Bei den gelesenen Namen, die möglicherweise "XA: Z:" enthalten, besteht eine Möglichkeit, Robustheit zu erzielen, darin, einen geeigneten Parser im SAM-Format wie pysam zu verwenden, wenn man in Python programmiert, anstatt grep. Ich weiß nicht, wie viel langsamer das sein würde.
Oder wie `` `samtools view```, das ich ausprobiert habe, aber nicht getan habe, was ich erwartet hatte. Ich kann mir vorstellen, dass samtools in Zukunft so etwas wie einen Tag-Filter implementieren könnte. Vielleicht hat es schon und ich benutze nur das falsche Sub-Tool.
Okay, es ist noch nicht in samtools verfügbar, aber die Idee, Tag-Filterung durchzuführen, wurde bereits in mindestens [einem Problem] angesprochen (https://github.com/samtools/samtools/issues/665#issuecomment-293926546).
Im schlimmsten Fall können Sie ein paar Zeilen Python schreiben oder wahrscheinlich Pierre's Jvarkit verwenden, um Dinge herauszufiltern.
#2
+3
Samir
2018-06-22 07:09:02 UTC
view on stackexchange narkive permalink

Für eine viel schnellere und stabilere $ sup> -Implementierung kann Sambamba dies mit dem folgenden Einzeiler handhaben:

  Sambamba-Ansicht -t 12 -h -f bam -F "Mapping_Quality > = 1 und nicht (nicht zugeordnet oder sekundäre Ausrichtung) und nicht ([XA]! = null oder [SA]! = null)" test.bwa-mem.bam -o test-uniq.bam  

Weitere Informationen zur weiteren Verwendung finden Sie unter manpage und sambamba synatx für Filter.

$: Neben grep Da die Suche langsam ist, werden möglicherweise falsche positive del> negative Treffer zurückgegeben (Antwort von @ gringer oben). Ich finde den awk -basierten Ansatz auch nicht zuverlässig, da Felder wie XA, SA usw. optionale Felder im SAM-Format sind und nicht positionell festgelegt sind Felder.

Inwiefern waren die Treffer falsch positiv?
Hmm .. Ich denke, es sollte falsch negativ sein, weil Sie diese Suchmuster mit `grep -v` ausgeschlossen haben, oder?
Die einzige Möglichkeit, wie grep zu falschen Ergebnissen führen kann, besteht darin, dass jemand die gelesenen Namen oder eine Lesegruppe so ändert, dass sie eine dieser Zeichenfolgen enthält, was ziemlich unwahrscheinlich erscheint.
stimme zu, es ist ziemlich unwahrscheinlich, aber dennoch ist der grep-basierte Ansatz von Natur aus langsam.
Der langsame Teil ist die Konvertierung von BAM nach SAM durch "samtools view". Das "grep" -Programm ist sehr schnell in dem, was es tut.


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