Zählvorkommen um eine genomische Region in einem Datenrahmen

Question

Ich muss Mutationen im Genom zählen, die an bestimmten Stellen oder Bereichen auftreten. Die Mutationen haben eine genomische Position (Chromosomen- und Basenpaar, z.B. Chr1, 10658324). Der Bereich bzw. der Spot ist definiert als 10000 Basenpaare aufwärts und abwärts (+ -) einer gegebenen Position im Genom. Sowohl Positionen von Mutationen als auch die Position von "Spots" sind in Datenrahmen gespeichert.

Beispiel:

set.seed(1) 

Chr <- 1 
Pos <- as.integer(runif(5000 , 0, 1e8)) 
mutations <- data.frame(Pos, Chr) 

Chr <- 1 
Pos <- as.integer(runif(50 , 0, 1e8)) 
spots <- data.frame(Pos, Chr) 

So ist die Frage ich frage ist: Wie viele Mutationen vorhanden sind + -10K Basenpaare um die gegebenen Positionen in "Spots". (z.B. wenn der Spot 100k ist, wäre der Bereich 90k-110k). Die echten Daten würden natürlich alle 24 Chromosomen enthalten, aber der Einfachheit halber können wir uns vorerst auf ein Chromosom konzentrieren. Die endgültigen Daten sollten den "Punkt" und die Anzahl der Mutationen in seiner Umgebung enthalten, idealerweise in einem Datenrahmen oder einer Matrix.

Vielen Dank im Voraus für Anregungen oder Hilfe!

---

Hier ist ein erster Versuch, aber ich bin ziemlich sicher, es gibt einen Weg, eleganter Weise, es zu tun.

w <- 10000 #setting range to 10k basepairs 
loop <- spots$Pos #creating vector of positions to loop through 
out <- data.frame(0,0) 
colnames(out) <- c("Pos", "Count") 

for (l in loop) { 
    temp <- nrow(filter(mutations, Pos>=l-w, Pos<=l+w)) 
    temp2 <- cbind(l,temp) 
    colnames(temp2) <- c("Pos", "Count") 
    out <- rbind(out, temp2) 
} 
out <- out[-1,] 

Answer 1

data.table foverlaps verwenden, dann Aggregat:

library(data.table) 
#set the flank 
myFlank <- 100000 

#convert to ranges with flank 
spotsRange <- data.table(
    chr = spots$Chr, 
    start = spots$Pos - myFlank, 
    end = spots$Pos + myFlank, 
    posSpot = spots$Pos, 
    key = c("chr", "start", "end")) 

#convert to ranges start end same as pos 
mutationsRange <- data.table(
    chr = mutations$Chr, 
    start = mutations$Pos, 
    end = mutations$Pos, 
    key = c("chr", "start", "end")) 

#merge by overlap 
res <- foverlaps(mutationsRange, spotsRange, nomatch = 0) 

#count mutations 
resCnt <- data.frame(table(res$posSpot)) 
colnames(resCnt) <- c("Pos", "MutationCount") 
merge(spots, resCnt, by = "Pos") 
#   Pos Chr MutationCount 
# 1 3439618 1   10 
# 2 3549952 1   15 
# 3 4375314 1   11 
# 4 7337370 1   13 
# ... 

Answer 2

Ich bin nicht mit Bett Manipulationen in R vertraut, also werde ich hier kann eine Antwort mit bedtools und jemand vorschlagen, versuchen, Granges oder andere R Bioinformatik Bibliothek zu konvertieren.

Im Wesentlichen haben Sie zwei Bett-Dateien, eine mit Ihren Flecken und andere mit Ihren Mutationen (ich nehme eine 1bp-Koordinate für jede in letzterem an). In diesem Fall würden Sie closestBed verwenden, um den nächstgelegenen Punkt und den Abstand in bp jeder Mutation zu ermitteln und dann diejenigen zu filtern, die 10 KB von den Punkten entfernt sind. Der Code in einer UNIX-Umgebung würde wie folgt aussehen:

# Assuming 4-column file structure (chr start end name) 
closestBed -d -a mutations.bed -b spots.bed | awk '$9 <= 10000 {print}' 

Wo Spalte 9 ($9) wird der Abstand in bp vom nächsten Ort und Stelle sein. Je nachdem, wie spezifisch Sie sein möchten, können Sie die Handbuchseite unter http://bedtools.readthedocs.io/en/latest/content/tools/closest.html überprüfen. Ich bin mir ziemlich sicher, dass es mindestens ein Betttools-ähnliches Paket in R gibt. Wenn die Funktionalität ähnlich ist, können Sie genau dieselbe Lösung anwenden.

Hoffe, dass hilft!