In der Genomik muss man oft die nächstgelegenen Varianten (z. B. SNPs, Indels) zu einer Reihe interessanter Genomorte finden. Das geht mit verschiedenen Bioinformatik-Werkzeugen wie bedtools. In diesem Beitrag zeige ich Schritt für Schritt, wie das funktioniert.
Voraussetzungen
Bevor es losgeht, sollten folgende Dateien vorhanden sein:
- Eine BED-Datei mit den interessierenden Orten. In diesem Beispiel nehmen wir
locations_of_interest.bed - Eine VCF-Datei mit den Varianten. In diesem Beispiel nehmen wir
FinalSetVariants_referenceGenome.vcf
Schritt 1: Die VCF-Datei sortieren
Das erste Problem war, dass die VCF-Datei nicht lexikografisch sortiert war. bedtools verlangt aber genau diese Sortierung. Die VCF-Datei lässt sich mit folgendem Befehl sortieren:
(grep '^#' FinalSetVariants_referenceGenome.vcf; grep -v '^#' FinalSetVariants_referenceGenome.vcf | sort -k1,1 -k2,2n) > sorted_FinalSetVariants_referenceGenome.vcf
Der Befehl trennt die Kopfzeilen (die mit # beginnen) von den Datenzeilen, sortiert die Datenzeilen und fügt Kopf und sortierte Daten in einer neuen Datei sorted_FinalSetVariants_referenceGenome.vcf zusammen.
Schritt 2: VCF in BED umwandeln und die nächsten Varianten finden
Als Nächstes suchen wir die nächstgelegenen Varianten zu unseren Orten. Standardmäßig gibt bedtools closest aber den ganzen VCF-Eintrag aus, was oft mehr ist als nötig. Um die Ausgabe zu begrenzen, wandeln wir die VCF-Datei unterwegs ins BED-Format um und vergeben als zusätzliches Merkmal den Markernamen als chr_bpLocation (so benennen wir unsere Marker). Mit der Option -d bekommen wir außerdem den Abstand zwischen dem Ort und der nächsten Variante. Hier der Befehl:
awk 'BEGIN {OFS="\t"} {if (!/^#/) {print $1,$2-1,$2,$4"/"$5,"+",$1"_"$2}}' sorted_FinalSetVariants_referenceGenome.vcf | bedtools closest -a locations_of_interest.bed -b stdin -d
Der Befehl liest mit awk die VCF-Daten, wandelt sie ins BED-Format um und schreibt das Ergebnis auf die Standardausgabe. Die Pipe (|) gibt diese Ausgabe direkt als -b-Datei an bedtools closest weiter. Das Schlüsselwort stdin sagt bedtools, dass es von der Standardeingabe lesen soll.
Fazit
Mit diesen zwei Schritten findet man effizient die nächstgelegenen Varianten zu einer Reihe interessanter Genomorte. Der Ansatz ist flexibel und lässt sich an andere Datensätze und Anforderungen anpassen.
