Frage:
Wie filtere ich Kreuzausrichtungen aus einer BED-Datei heraus?
SmallChess
2017-05-19 10:49:47 UTC
view on stackexchange narkive permalink

Ich habe eine BAM-Datei:

  @SQ SN: chr1 LN: 248956422 @ SQ SN: chrx LN: 248956423ST-E00110: 348: HGVKKALXX: 1: 1201: 5822: 48670 323 chr1 9999 0 67H66M16H CHRX 1000 0 GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ NM: i: 0 MD: Z: 66 AS: i: 66 XS: i: 65 SA: Z: chr5,18606834, -, 73S76M, 34,0; RG: Z: g1  

Es gibt einen Lesevorgang, der auf chr1 ausgerichtet ist, und einen Partner, der auf chrx ausgerichtet ist.

Ich habe eine BETT-Datei:

  chr1 0 100000 TestOnly  

Ich möchte alles herausfiltern, was außerhalb meines BETTES liegt Region, die Kreuzausrichtungen enthält. In meinem Beispiel ist mein Lesevorgang zwar auf chr1 ausgerichtet, aber es ist kein Partner. Ich möchte nicht, dass dies gelesen wird.

Wenn ich dies tue:

samtools-Ansicht -L test.bed test.bam

Der Befehl gibt mir den Lesevorgang, da keine Kreuzausrichtungen überprüft werden.

Meine Lösung:

samtools view -L test.bed test.bam | grep -v chrx

aber das ist sehr langsam und ungeschickt. In meiner Produktionspipeline müsste ich Folgendes tun:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

F: Gibt es eine bessere Lösung?

Einer antworten:
#1
+6
terdon
2017-05-19 22:44:29 UTC
view on stackexchange narkive permalink

Gemäß der SAM-Spezifikation lautet das dritte Feld einer SAM-Zeile ( RNAME ):

RNAME: Referenzsequenz NAME der Ausrichtung. Wenn @ SQL-Headerzeilen vorhanden sind, muss RNAME (falls nicht "*") in einem der SQ-SN-Tags vorhanden sein. Ein nicht zugeordnetes Segment ohne Koordinate hat in diesem Feld ein "*". Ein nicht zugeordnetes Segment kann jedoch auch eine gewöhnliche Koordinate haben, so dass es nach dem Sortieren an einer gewünschten Position platziert werden kann. Wenn RNAME '*' ist, können keine Annahmen über POS und CIGAR getroffen werden.

Und das 7. Feld lautet (Hervorhebung von mir, es fehlt "bis"):

RNEXT: Name der Referenzsequenz der primären Ausrichtung des in der Vorlage gelesenen NEXT. Beim letzten Lesevorgang ist der nächste Lesevorgang der erste Lesevorgang in der Vorlage. Wenn @ SQ-Kopfzeilen vorhanden sind, muss RNEXT (wenn nicht "*" oder "=") in einem der SQ-SN-Tags vorhanden sein. Dieses Feld wird als "*" festgelegt, wenn die Informationen nicht verfügbar sind, und als "=", wenn RNEXT mit RNAME identisch ist. Wenn nicht "=" und der nächste Lesevorgang in der Vorlage eine primäre Zuordnung aufweist (siehe auch Bit 0x100 in FLAG), ist dieses Feld mit RNAME in der primären Zeile des nächsten Lesevorgangs identisch. Wenn RNEXT '*' ist, können für PNEXT und Bit 0x20

keine Annahmen getroffen werden. Sie möchten also die Zeilen entfernen, deren 7. Feld nicht = ist und nur für den Fall, dass die Zeilen, deren 7. Feld nicht = und ist, nicht mit dem 3. Feld identisch sind. Sie können daher Folgendes verwenden:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" || $ 3 == $ 7  

Und um erneut als BAM-Datei zu speichern:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" && $ 3 == $ 7 | samtolls view -b > fixed.bam  

In einem separaten Hinweis ist es sehr selten erforderlich, mehrere solche grep-Befehle zu verketten. Sie können sie einfach mit \ | (oder | mit den Optionen -E oder -P ) trennen. So etwas wie:

  samtools view -L test.bed test.bam | grep -v 'chrx \ | chr2 \ | chr10 \ | chrN'  

Oder

  samtools view -L test.bed test.bam | grep -Ev 'chrx | chr2 | chr10 | chrN'  
Wenn Sie dies so tun, fehlt in der Datei "fixed.bam" der Header, was meiner Erfahrung nach viele Probleme verursacht. Ich empfehle, den Header immer wieder hinzuzufügen. entweder durch Angabe von "-h" beim Lesen der ursprünglichen BAM oder durch separates Hinzufügen: "(samtools view -H infile.bam; samtools view ...)> samtools view -b> outfile.bam".


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