Frage:
Untermenge kleinerer BAM, um mehrere tausend Zeilen von mehreren Chromosomen zu enthalten
EB2127
2018-02-17 18:00:18 UTC
view on stackexchange narkive permalink

Es gibt viele Fälle, in denen ich eine BAM unterteilen möchte, um eine kleine Datei zu erstellen, mit der gearbeitet werden kann (z. B. algorithmisches Testen, Debuggen usw.).

Normalerweise mache ich Folgendes Teilmenge der BAM file.bam und Behalte den Header

  samtools view -H file.bam > header.samsamtools view file.bam | Kopf -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam  

In diesem Fall möchte ich 5000 Zeilen in Chromosom 1 und 5000 in Chromosom 2.

Ich könnte zuerst Versuchen Sie, nach einzelnen Chromosomen zu greifen, und kombinieren Sie dann die beiden SAMs

, z Hier ist die vollständige BAM mit grepped chr1 und (falschem, aber vollständigem) Header

  samtools view -H file.bam > header.samsamtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam  

aber dann habe ich zwei Probleme:

(1) Möglicherweise greife ich die Ausrichtungen zu Chromsome 2- nicht. - Es können BAM-Zeilen vorhanden sein, die 'chr2' enthalten, aber keine Alignments sind.

(2) Ich denke, man muss den Header manuell bearbeiten. Daran führt wahrscheinlich kein Weg vorbei.

Gibt es einen einfachen Weg, Bioinformatik SO?

Fünf antworten:
Devon Ryan
2018-02-17 18:29:12 UTC
view on stackexchange narkive permalink

Wenn Sie nicht zu sehr auf genaue Zahlen wie 5000 Lesevorgänge angewiesen sind, können Sie dies mit einem einzigen samtools-Befehl tun:

  samtools-Ansicht -bo subset.bam -s 123.4-Ausrichtungen. bam chr1 chr2  

Damit werden 40% (der .4 -Teil) der Lesevorgänge ausgewählt ( 123 ist ein Startwert) bequem für die Reproduzierbarkeit). Der praktische Teil davon ist, dass die Partner gepaart bleiben, wenn Sie Lesevorgänge am gepaarten Ende haben. Ändern Sie für 5000 Lesevorgänge pro Chromosom einfach den Teil .4 in eine ausreichend kleine Zahl.

Im Allgemeinen müssen Sie den Header nicht wirklich unterteilen. Einige Tools erzielen in diesem Fall eine etwas bessere Leistung, erzielen jedoch im Allgemeinen unabhängig davon dieselben Ergebnisse.

Karel Brinda
2018-02-27 04:30:15 UTC
view on stackexchange narkive permalink

Sie können SAMsift verwenden:

  samsift \ -i file.bam \ -0 'c = {"chr1": 5000, "chr2": 5000 } '\ -f' c [RNAME] >0 '\ -c' c [RNAME] - = 1 '\ -m nonstop-remove  

Erläuterung:

  • -i file.bam - Eingabedatei
  • -0 'c = {"chr1": 5000, "chr2": 5000}' - Initialisierung (Erstellen eines Countdown-Wörterbuchs für die interessierenden Chromosomen)
  • -f 'c [RNAME] >0' - Filterkriterium (ist der Zähler für das aktuelle Chromosom immer noch> 0?)
  • -c 'c [RNAME] - = 1' - Code, der den Zähler des aktuellen Chromosoms dekrementiert (5000 -> 4999 -> ... - > 1 -> 0 -> -1 -> ...)
  • -m nonstop-remove - Entfernen Sie Zeilen, die Python-Fehler verursachen, und hören Sie nicht auf (in diesem Fall, Ein Fehler kann durch den Zugriff auf einen nicht vorhandenen Zähler für ein anderes Chromosom verursacht werden, z. B. für chr3)

Weitere Informationen finden Sie in der SAMsift-Readme.

Wenn Sie eine Erklärung hinzufügen könnten, was diese Optionen sind und warum sie notwendig sind. Zum Beispiel verstehe ich nicht, warum das Argument für "f" "c [RNAME]" ist und nicht andere Dinge
Ich habe gerade die Erklärung hinzugefügt.
Ian Sudbery
2018-02-17 21:20:09 UTC
view on stackexchange narkive permalink

Ich würde Devon Ryans Antwort generell empfehlen. Wenn Sie jedoch die gleiche Nummer von jedem Chr haben möchten, können Sie den folgenden Python / Pysam-Code verwenden (dies gibt ungefähr 5000 von jedem Chr aus):

  aus pysam import AlignmentFilefrom random import randomnreads = 5000infile = AlignmentFile (" mybam.bam ") outfile = AlignmentFile (" outbam.bam "," wb ", template = infile) für chr in infile.reference_names: read_in_chr = infile.count (chr) frac = float (nreads) / read_in_chr count = 0 # - Ersetzen Sie diese Schleife für genau 5000 erste Lesevorgänge - # für das Einlesen von infile.fetch (chr, multiple_iterators = True): if random () < frac: count + = 1 outfile.write (read) print ("ouputted% i liest aus% s"% (count, chr)) outfile.close ()  

Dies verwendet alle Chromosomen. Wenn Sie nur chr1 und chr2 verwenden möchten, ersetzen Sie infile.reference_names durch ["chr1","chr2"‹.

Wenn Sie möchten genau 5000 Lesevorgänge von jedem chr, aber egal ob sie die ersten sind oder nicht, dann können Sie die innere for-Schleife durch Folgendes ersetzen:

  zum Einlesen von infile.fetch (chr, multiple_iterators = True): count + = 1, wenn count > = nreads: break outfile.write (read)  

, wenn Sie die Partner dieser haben möchten Wenn Sie auch lesen, können Sie nach outfile.write (read) Folgendes hinzufügen:

  mate_read = infile.mate ( read) outfile.write (mate_read)  

Beachten Sie, dass dies langsam ist.

Um die Arbeit zu beschleunigen, kann "infile.count (chr)" durch das Parsen von "idxstats ()" ersetzt werden (es sei denn, man verwendet CRAM-Dateien).
Vielen Dank, ich hatte angenommen, dass es eine einfache Möglichkeit gibt, dies aus dem Index abzurufen.
"Um die Dinge zu beschleunigen, kann infile.count (chr) durch das Parsen von idxstats () ersetzt werden." Ich folge nicht ganz. Wie würde das oben genannte aussehen?
Pierre
2018-02-19 19:49:23 UTC
view on stackexchange narkive permalink

using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

  $ java -jar dist /samjdk.jar --body -e \ 'Map<String, Integer> c = new HashMap<> (); public Object apply (SAMRecord r) {int n = c.getOrDefault (r.getContig (), 0); if (n> = 5000) return false; c.put (r.getContig (), n + 1); return r;} '\ input.bam  
Es wäre schön, wenn die Leute, die ihre Antworten ablehnen, erklären würden, warum sie das tun ...
@Devon (ich habe den Kommentar bis jetzt gesehen) Ich habe abgelehnt, weil ich finde, dass eine lange Codezeile, ohne zu erklären, was sie intern tut, normalerweise schwerer zu verstehen ist. Der Code, wie er sich verbessern würde (IMHO), ist in mehreren Zeilen und mit Kommentaren gebrochen . Es fällt mir schwer, dass ich Java nicht kenne, um zu verstehen, was passiert
Konrad Rudolph
2018-02-20 21:08:42 UTC
view on stackexchange narkive permalink

Wenn Sie nur eine abgeschnittene BAM-Datei mit einem Header erstellen möchten, können Sie Ihren ursprünglichen Code erheblich vereinfachen:

  (samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam  

Dieser Befehl vermeidet die Zwischendatei und einige kostenlose Aufrufe von Zwischenbefehlen auf Kosten einer Unterschale (aufgerufen durch (…) ). .

Wie oben, jedoch mit Formatierung:

  (samtools-Ansicht -H input.bam samtools-Ansicht input.bam | head -5000 # (*)) \ | samtools -bo output.bam  

… dies kann natürlich erweitert werden, um durch mehrere Chromosomen zu filtern, indem die mit (*) oben markierte Zeile durch eine oder mehrere ersetzt wird Zeilen, die nach Chromosomennamen untergeordnet sind ( samtools view input.bam chr x , grep ist nicht erforderlich, wenn Sie die ursprüngliche BAM-Datei indiziert haben!) .

Ist es `samtools -bo`?


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