Frage:
Wie kann ich das Aufrufen / Korrigieren von INDEL für BAM-Dateien beschleunigen?
gringer
2017-06-05 04:25:12 UTC
view on stackexchange narkive permalink

Der Befehl samtools mpileup verfügt über eine nette Funktion, mit der Zuordnungsfehler korrigiert werden können, die mit einer Fehlausrichtung von INDELs verbunden sind. Standardmäßig funktioniert der Befehl mpileup nicht für Lesevorgänge, bei denen das Referenzgenom mehr als 250-fach abgedeckt ist. Obwohl dieses Limit erhöht werden kann, führt eine sehr hohe Abdeckung dazu, dass das mpileup-Programm zum Stillstand kommt. Es wäre also schön zu wissen, ob es eine einfache Möglichkeit gibt, dies zu beschleunigen.

Um ein bisschen mehr hinzuzufügen Im Kontext habe ich dies mit mitochondrialen Genom-Reads gemacht, die sowohl aus der Illumina-Gesamtgenomsequenzierung (Abdeckung ~ 1000X) als auch aus der gezielten Amplikonsequenzierung auf dem IonTorrent (Abdeckung bis ~ 4000X) extrahiert wurden.

Ich sehe, dass @rightskewed die Downsampling-Fähigkeit von samtools mit samtools view -s <float> erwähnt hat (siehe hier), was als Lösung dafür zu gelten scheint wenn vor dem mpileup-Vorgang verwendet.

Einer antworten:
#1
+4
gringer
2017-06-05 04:34:22 UTC
view on stackexchange narkive permalink

Als ich vor ein paar Jahren dieses Problem hatte, war mir das samtools-Subsampling nicht bekannt. Daher schrieb ich meine eigene digitale Normalisierungsmethode, um mit zugeordneten Lesevorgängen umzugehen. Diese Methode reduziert die Genomabdeckung, behält jedoch Lesevorgänge bei geringer Abdeckung bei.

Da ich mit IonTorrent-Lesevorgängen (mit variabler Länge) gearbeitet habe, kam mir die Idee, den längsten Lesevorgang auszuwählen, dem zugeordnet wurde jeder Ort im Genom (vorausgesetzt, ein solcher Lesevorgang existiert). Dies bedeutete, dass die sehr variable Abdeckung für verschiedene Proben (manchmal so niedrig wie 200X, manchmal so hoch wie 4000X) auf eine viel konsistentere Abdeckung von ungefähr 100-200X abgeflacht wurde. Hier ist der Kern des Perl-Codes, den ich geschrieben habe:

  if (($ F [2] ne $ seqName) || ($ F [3]! = $ Pos) || (Länge ($ bestSeq) < = Länge ($ F [9])) {if (Länge ($ bestSeq) == Länge ($ F [9])) {## Reservoir-Probenahme mit einer Reservoirgröße von 1 ## Siehe https : //en.wikipedia.org/wiki/Reservoir_sampling ## * mit Wahrscheinlichkeit 1 / i, behalte das neue Element anstelle des aktuellen Elements $ foundCount ++; if (! rand ($ SeenCount)) {## d. h. wenn rand ($ SeenCount) == 0 ist, fahren Sie mit dem nächsten Ersetzen fort; }} else {$ foundCount = 1; } if (($ F [2] ne $ seqName) || ($ F [3]! = $ pos)) {if ($ output eq "fastq") {printSeq ($ bestID, $ bestSeq, $ bestQual); } elsif ($ output eq "sam") {print ($ bestLine); }} $ seqName = $ F [2]; $ pos = $ F [3]; $ bestLine = $ line; $ bestID = $ F [0]; $ bestFlags = $ F [1]; $ bestSeq = $ F [9]; $ bestQual = $ F [10];  

Der vollständige Code (der als Filter für unkomprimierte SAM-Dateien fungiert) finden Sie hier.

Zur Verdeutlichung: Dies ist Perl-Code?
Ja, Perl-Code. Ich kann es überall an den Dollarsymbolen erkennen, aber ich denke, es ist für Leute, die Perl nicht kennen, nicht offensichtlich.
Nun, ich spreche Perl, aber es ist kaum die einzige Sprache mit Dollar vor Identifikatoren. ;-);
Tatsächlich. Meine Gedanken waren leer, aber meine Frau erinnerte mich an Bash / Sh.


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