教程:
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添加辅助线
如何在已有图形上面添加一条直线?
使用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,颜色为红色
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"))