④改进版 valr实例计算ChIP-seq中Upstream +

数据导入：这里使用valr包自带的数据

library(valr)
bedfile <- valr_example("genes.hg19.chr22.bed.gz")  ## 基因或者目标区域坐标信息
genomefile <- valr_example("hg19.chrom.sizes.gz")  ## 染色体长度信息
bgfile  <- valr_example("hela.h3k4.chip.bg.gz")  ## 信号值文件

genes <- read_bed(bedfile, n_fields = 6)
genome <- read_genome(genomefile)
y <- read_bedgraph(bgfile)

导入前的数据格式为

image.png

写成函数, 计算区域内的信号值

get_right_Gene_bins_cov <- function(bedfile, num_win, y){
## bedfile 即我们的基因区域文件
## num_win 表示我们需要划分成多少个bin
## y  表示我们需要的信号值
  F_Genes <- bedfile %>%
    filter(strand == "+") %>%
    bed_makewindows(num_win = num_win)  
## 筛选正向的坐标，如果为正向，可以直接正常分
  
  R_Genes <- bedfile %>%
    filter(strand == "-") %>%
    bed_makewindows(num_win = num_win, reverse = T)
## 如果是反向的，那么最后分组的顺序倒置过来
## 操作见之前的valr说明

  new_Genes <- rbind(F_Genes, R_Genes) 
## 合并得到的正向和反向的文件
  GeneBody <- bed_map(new_Genes, y, win_sum = sum(value, na.rm = TRUE)) %>%
## bed_map 即根据所提供的信号值文件y, 计算我们目标分好的每一个bin中的信号值的总数
    group_by(.win_id) %>% 
## 按照我们划分bin时候，bed_makewindows生成的20个bin的编号来分组
    summarize(win_mean = mean(win_sum, na.rm = TRUE))
## 计算每一个bin中对应的所有区域数据的信号值的平均值
}

通过函数计算genebody区的信号值，将其划分为20个bin，最后的提取信号值列

GeneBody <- get_right_Gene_bins_cov(genes, 20, y)[,2]  

计算Promoter区每一个bin中的信号

region_size <- 2000
# 50 bp windows
win_size <- 100
bin_number <- region_size *2 /win_size
## 即 Promoter区域加downstream区域总共40个bin

#### 获取Promoter的信息
F_promoter <- genes %>%
  filter(strand == "+") %>%
  bed_flank(genome, left = region_size) %>%
  get_right_Gene_bins_cov(20, y)

R_promoter <- genes %>%
  filter(strand == "-") %>%
  bed_flank(genome, right = region_size) %>%
  get_right_Gene_bins_cov(20, y)

## 再将正负信息合并，取平均值
Promoter <- data.frame(win_mean = apply(merge(F_promoter, R_promoter, by = ".win_id")[,2:3], 1, mean))

计算downstream区每一个bin中的信号

#### 获取down的信息
F_down <- genes %>%
  filter(strand == "+") %>%
  bed_flank(genome, right = region_size) %>%
  get_right_Gene_bins_cov(20, y)

R_down <- genes %>%
  filter(strand == "-") %>%
  bed_flank(genome, left = region_size) %>%
  get_right_Gene_bins_cov(20, y)

## 再将正负信息合并，取平均值
Down <- data.frame(win_mean = apply(merge(F_down, R_down, by = ".win_id")[,2:3], 1, mean))

合并Promoter、GeneBody、Downstream三个区域的信息，得到一个20*3 = 60行的矩阵

bind_cov <- rbind(Promoter, GeneBody, Down)
end_cov <- cbind(.win_id = seq_len(nrow(bind_cov)), bind_cov)
# 即前20行为Promoter区，中间20行为GeneBody区，最后20行为Downstream区

得到了矩阵就很方便了，使用ggplot2画图

library(ggplot2)
ggplot(end_cov, aes(x = .win_id, y = win_mean)) +
  geom_point(size = 0.25) + geom_line(size = 1)  +
  scale_x_continuous(labels = c("-2Kb","TSS", "TES", "2Kb"), breaks = seq(0,bin_number*3/2, by = bin_number/2)) +
  theme_bw() + xlab("") + ylab("Signal") +
  theme(legend.key=element_rect(linetype='dashed',color="white"),
        axis.text.y = element_text(size=20),axis.text.x = element_text(size=20),
        legend.title = element_text(size=15),legend.text = element_text(size=15),
        legend.key.height=unit(1.2,'cm')) 

image.png

将信号值取log2进行可视化

ggplot(end_cov, aes(x = .win_id, y = log2(win_mean))) +
  geom_point(size = 0.25) + geom_line(size = 1)  +
  scale_x_continuous(labels = c("-2Kb","TSS", "TES", "2Kb"), breaks = seq(0,bin_number*3/2, by = bin_number/2)) +
  theme_bw() + xlab("") + ylab("Signal") +
  theme(legend.key=element_rect(linetype='dashed',color="white"),
        axis.text.y = element_text(size=20),axis.text.x = element_text(size=20),
        legend.title = element_text(size=15),legend.text = element_text(size=15),
        legend.key.height=unit(1.2,'cm')) 

image.png

最后补充关于bed_makewindows函数,举例

> x <- trbl_interval(~chrom, ~start, ~end, ~name, ~score, ~strand,
                   "chr1", 100,    200,'A','.','+',
                   "chr1", 300, 400, 'B', '.', '-')
> bed_makewindows(x, num_win = 5)
# A tibble: 10 x 7
   chrom start   end name  score strand .win_id
   <chr> <int> <int> <chr> <chr> <chr>    <int>
 1 chr1    100   120 A     .     +            1
 2 chr1    120   140 A     .     +            2
 3 chr1    140   160 A     .     +            3
 4 chr1    160   180 A     .     +            4
 5 chr1    180   200 A     .     +            5
 6 chr1    300   320 B     .     -            1
 7 chr1    320   340 B     .     -            2
 8 chr1    340   360 B     .     -            3
 9 chr1    360   380 B     .     -            4
10 chr1    380   400 B     .     -            5
> bed_makewindows(x, num_win = 5, reverse = T)
# A tibble: 10 x 7
   chrom start   end name  score strand .win_id
   <chr> <int> <int> <chr> <chr> <chr>    <int>
 1 chr1    100   120 A     .     +            5
 2 chr1    120   140 A     .     +            4
 3 chr1    140   160 A     .     +            3
 4 chr1    160   180 A     .     +            2
 5 chr1    180   200 A     .     +            1
 6 chr1    300   320 B     .     -            5
 7 chr1    320   340 B     .     -            4
 8 chr1    340   360 B     .     -            3
 9 chr1    360   380 B     .     -            2
10 chr1    380   400 B     .     -            1

可以看出来，在不指定reverse参数为True时候，只有正方向的bin的划分格式为正确的，反方向的bin划分的时候应该为加了reverse参数为True时候这种划分方式，即应该是反的。所以我们区分正负方向划分bin的时候只需先将正负方向各自提取出来，然后正方向按照正常的方式来划分bin，反方向按照reverse反方向划分bin，然后在合并两个文件的bin的信息，同一个bin取均值即可。

同理，我们在得到启动子和下游2kb区间时候利用bed_flank函数的right 和 left参数，划分好正负反向然后合并，即可。

image.png