Rolling_average_plots

教程：

bedtools的使用技巧（持续更新）
bedtools使用指南
bedtools统计窗口内平均覆盖深度

• 目的：
  画一个类似这样的图
  

思路：

1.划分参考基因组，将参考基因组划分为大小一样的n个窗口。比如划分hg19，每个窗口大小10kb/0.5Mbp。将这些窗口作为后续计算reads数的区间。
2.利用bedtools coverage 计算每个区间的reads数，并根据RRPM进行normalization。
3.利用metaplot/python的时间序列分析（pandas函数）进行作图。

1.划分bins

利用bedtools makewindows 划分bins

 bedtools makewindows -g /Data/wu_lab/reference/UCSC_hg19/hg19.chrom.sizes.1 -w 500000 > windows.bed > windows.bed

-g:参考基因组 格式为：染色体号，染色体长度
如：chr8 xxxxxxx
-w：划分为大小为0.5M的窗口

去掉hg19的基因黑名单

bedtools intersect -wa -wb -a windows.bed -b dukeExcludeRegions.bed -v > hg19_no_blacklist.bed

2.计算reads （统计窗口内的平均覆盖）

#!/bin/bash

for i in DIPG17-combine-H3K27me3 DIPG17-combine-H3K27me3_dm DIPG17-control-H3K27me3 DIPG17-control-H3K27me3_dm DIPG17-EC8042-H3K27me3 DIPG17-EC8042-H3K27me3_dm DIPG17-Vorinostat-H3K27me3 DIPG17-Vorinostat-H3K27me3_dm

do 
bedtools coverage -a windows.bed -b ${i}.bam.sort > ${i}.xls

done

-a a指定interval文件，即你想看的染色体区间
-b指定的是你比对的结果（bam）或bed等文件

得到以下文件格式：
本次数据做了normalization，因此画图前需要算出每个样本的scale factor
将每个excel整理成以下格式
提取control以及每个样本的reads 。然后用如EC8042_reads/control_reads，×scale factor（每组数值不同）。并取log2。整理成下表格式。

3.画图

setwd("F:/chip-seq/kongyu/2022_4_5_dipg_drug/bin_p2/")
library(ggplot2)
data_EC8042_bin <- read.csv("log2FC_EC8042_vs_control_bins0.5M_chr8.csv")
head(data_EC8042_bin)
x=data_EC8042_bin$bins
y=data_EC8042_bin$log2FC_EC8042_control
head(x)
plot(x,y,type = "l",col="blue",ylim=c(-2,2),xlab = "chr8",ylab="log2FC_EC8042_control", main = "EC8042 vs control H3K27me3",xaxt="n")>abline(h=0,col="black") 

结果图

plot添加辅助线

R语言在已有图形上面添加一条直线

如何在已有图形上面添加一条直线？
使用abline()函数。示例代码如下。
>x <- c(1.0, 2.0, 3.0, 4.0, 5.0)
>y <- c(1.0, 1.9, 3.1, 4.0, 4.9)
>plot(x,y)>abline(v=3,lwd=4,col="blue")#添加一条垂直直线x=3，线宽为4，颜色蓝色
>abline(h=3,lwd=4,col="blue")#添加一条水平直线y=3，线宽为4，颜色蓝色
>abline(lm(y~x), lwd=4, col="red")#添加一条一元线性回归拟合直线，线宽为4，颜色为红色

R技巧[39]：R语言绘图修改坐标轴刻度

plot(yt,type="l",xaxt="n"); #xaxt="n"表示去掉原来刻度

axis(1,c(1,43,80,122),labels=c("2018/11","2019/01","2019/03","2019/05"))