How to create genomic windows in R?
I have methylation data which looks like this for example.
chr start end meth
chr7 31441 31441 0.16542433
chr7 31467 31467 0.93508003
chr7 50060 50060 0.38091076
chr7 50097 50097 0.31270269
chr7 50147 50147 0.39961491
chr7 50158 50158 0.01449239
chr7 50164 50164 0.76305873
chr7 50355 50355 0.56224390
chr7 75862 75862 0.04551076
chr7 79874 79874 0.57058403
I would like to create 1000bp window and take median of meth. It looks like following
window chr median start stop
0 chr7 NA 0 10000
1 chr7 NA 10000 20000
2 chr7 NA 20000 30000
3 chr7 0.5502522 30000 40000
4 chr7 NA 40000 50000
5 chr7 0.3902628 50000 60000
6 chr7 NA 60000 70000
7 chr7 0.3080474 70000 80000
I was able to create such table in R, but kind of incomplete. R Code:
df <- data.frame(matrix(ncol = 4, nrow = 10))
colnames(df) <- c("chr","start","end","meth")
df$chr <- c(rep("chr7",10))
df$start <- c(31441,31467,50060,50097,50147,50158,50164,50355,75862,79874)
df$end <- df$start
df$meth <- runif(10, min = 0, max = 1)
bin <- 10000
df2 <- df %>%
mutate(window = .$start %/% bin) %>%
group_by(window,chr) %>%
summarise(median = median(meth)) %>%
mutate(start = window*bin, stop=(window+1)*bin)
Output:
Source: local data frame [3 x 5]
Groups: window [3]
window chr median start stop
<dbl> <chr> <dbl> <dbl> <dbl>
1 3 chr7 0.5502522 30000 40000
2 5 chr7 0.3902628 50000 60000
3 7 chr7 0.3080474 70000 80000
• 3,664 views
•
link
updated
by
andrew.j.skelton73
•
written
by
Tori
•
written
by
Tori
1 answer
As far as I can see you've done the hard part, and you just need to fill in the gaps with the empty windows. See the below solution:
library(dplyr)
df <- data.frame(matrix(ncol = 4, nrow = 20))
colnames(df) <- c("chr","start","end","meth")
df$chr <- c(rep("chr7",10),rep("chr10",10))
df$start <- c(31441,31467,50060,50097,50147,50158,50164,50355,75862,79874,
31446,31467,50080,50099,50148,50160,50163,50352,75867,79879)
df$end <- df$start
df$meth <- runif(20, min = 0, max = 1)
bin <- 10000
df2 <- df %>%
mutate(window = start %/% bin) %>%
group_by(window,chr) %>%
summarise(median = median(meth)) %>%
mutate(start = (window*bin)+1, stop=(window+1)*bin)
df_out <- c()
for(i in unique(df2$chr)) {
df2_tmp <- df2 %>% filter(chr == i)
missing_win <- setdiff(0:max(df2_tmp$window),unique(df2_tmp$window))
df_tmp <- data.frame(window = missing_win) %>%
mutate(chr = i,
median = NA,
start = (window*bin)+1,
stop = (window+1)*bin) %>%
bind_rows(df2_tmp) %>%
arrange(chr,window)
df_out <- bind_rows(df_out,df_tmp)
}
df_out %>% arrange(desc(chr), window)
Results:
window chr median start stop
1 0 chr7 NA 1 10000
2 1 chr7 NA 10001 20000
3 2 chr7 NA 20001 30000
4 3 chr7 0.2974288 30001 40000
5 4 chr7 NA 40001 50000
6 5 chr7 0.5453479 50001 60000
7 6 chr7 NA 60001 70000
8 7 chr7 0.7934429 70001 80000
9 0 chr10 NA 1 10000
10 1 chr10 NA 10001 20000
11 2 chr10 NA 20001 30000
12 3 chr10 0.6825741 30001 40000
13 4 chr10 NA 40001 50000
14 5 chr10 0.5777722 50001 60000
15 6 chr10 NA 60001 70000
16 7 chr10 0.4252692 70001 80000
https://github.com/AndrewSkelton/Biostars-Answers/blob/master/p238223.Rmd
• 0 views
•
link
written
by
andrew.j.skelton73
Not the most elegant, but see the updated solution.
• 0 views
•
link
written
by
andrew.j.skelton73
Log in to answer this question.