ATAC-seq 片段长度周期性分布图

参考：

• 给你bam文件，你会画插入片段长度分布图吗？

• 糗世界-ATAC-seq

• 「文献02」esATAC: 一个ATAC-seq分析的R包

• 【r<-ggplot2】修改x和y轴刻度

• cowplot 小图嵌入大图

• 测试数据链接：CTCF_fragment.length.txt
  链接: https://pan.baidu.com/s/1WAdTqL8VUw84QDkc0qBisA 提取码: enyx

缘由：

经常看到ATAC 有一个片段长度周期性分布图，但是一直没注意右上角小图有什么含义？感觉是对fragment 统计的count 进行了log10，再画图
刚好最近跑了Tn5 酶切数据，就尝试了一下。

image.png

实践：

• Step1 ：获取insert fragment 长度数据

samtools  view input.bam | \
awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)} {print $1"\t"abs($9)}' | \
sort | uniq | cut -f2 > fragment.length.txt

• Step2 ： 画图

library(tidyverse)
data <- read.table("E:\\单细胞\\HZAU_scChIP\\测试CUT&Tag\\Fig1_insert_length\\CTCF_fragment.length.txt")

# 去除含零行
data <- tbl_df(data) %>% filter(V1!=0)
# 设置插入片段长度的阈值，过滤掉太长的片段
length_cutoff <- 1200
fragment <- data$V1[data$V1 <= length_cutoff]


######
##Part1：基础语法画图
######
# 利用直方图统计频数分布，设置柱子个数
breaks_num <- 500
res <- hist(fragment, breaks = breaks_num, plot = FALSE)


# 添加坐标原点
plot(x = c(0, res$breaks),
     y = c(0, 0, res$counts) / 10^2,
     type = "l", col = "red",full="red",
     xlab = "Fragment length(bp)",
     ylab = expression(Normalized ~ read ~ density ~ 10^2),
     main = "Sample Fragment sizes")

image.png

######
##Part2：ggplot2 画图及其拼图
######
## 不同数据分布
DATA <- data.frame(x1 = c(0, res$breaks),y1=c(0, 0, res$counts) / 10^2)
p1 <- ggplot(DATA,aes(x =x1,y = y1 ))+
  geom_line(col="red")+
  xlab("Fragment length(bp)")+
  ylab(expression(Normalized ~ read ~ count ~ 10^2))+
  ggtitle("Sample Fragment sizes")+
  theme_classic()
  
## 画小图
DATA2 <- data.frame(x1 = c(0, res$breaks),y1=log10(c(0, 0, res$counts) / 10^2)+1)
p2 <- ggplot(DATA2,aes(x =x1,y = y1 ))+
  geom_line(col="red")+
  xlab("Fragment length(bp)")+
  ylab(expression(Normalized ~ read ~ count ~ (log)))+
  ggtitle("Sample Fragment sizes")+
  theme_classic()


## 小图插入右上角
library(cowplot)
ggdraw() +
  draw_plot(p1, 0, 0, 1, 1) +
  draw_plot(p2, 0.5, 0.52, 0.5, 0.4) +
  draw_plot_label(c("A", "B"), c(0, 0.5), c(1, 0.92), size = 15)

image.png

######
##Part3：其他（密度图）
######
P3 <- ggplot(as.data.frame(fragment))+
  geom_density(aes(x=fragment),bw=3,color = "red")+
  xlim(0,NA)+
  scale_y_continuous(breaks = c(seq(0,0.004,0.001)),
                     labels=c(seq(0,4)),
                     name = expression(Normalized ~ read ~ density ~ 10^-3))+
  theme_classic()

修改成密度图

思考：

通过对insert fragment 取对数操作，可以知道右上角小图画法。为什么需要怎么操作呢？可能为了展示出密度图的波动性，及片段长度分布周期性。