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

欢迎关注”生信修炼手册”!web

对于ATAC文库而言，其插入片断的长度分布有着很是典型的规律，示意以下数据库

每200bp会存在一个峰，这个周期性波动反应的是核小体的个数。在ATAC_seq的数据分析中，会对插入片断长度分布进行可视化，观察其是否符合这样的周期性规律，必定程度能够反映文库构建的质量，那么如何在作这样一张分布图呢?express

比对以后咱们会获得bam文件，画图所需的插入片断长度就须要从bam文件中提取，须要注意，这里的插入片断是文库中adapter之间的插入片断，即fragment，须要和insert size区别开来。bash

对于单端测序而言，只有bam文件是不够的，须要借助工具来预测fragment length, 这里就不展开了。对于双端测序而言，事情就变得简单了。bam文件的第9列直接存储了fragment length的信息，直接提取以后就能够用来画图，提取的代码以下微信

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

bam文件中每一行以reads为单位，这里去重是为了不来自同一个fragemnts的reads重复统计。提取好以后，用R画图就能够了，R代码以下app

data <- read.table("fragment.length.txt", header = F)
# 设置插入片断长度的阈值，过滤掉太长的片断
length_cutoff <- 1200
fragment <- data$V1[data$V1 <= length_cutoff]
# 利用直方图统计频数分布，设置柱子个数
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",
 xlab = "Fragment length(bp)",
 ylab = expression(Normalized ~ read ~ density ~ 10^2),
 main = "Sample Fragment sizes")

输出结果示意以下less

这种是最简单的方式，除此以外，还有picard的CollectInsertSizeMetrics， bedtools的bamPEFragmentSize也均可以计算插入片断长度，可是在我看来，仍是本文介绍的方式最为简单直接。编辑器

·end·工具

—若是喜欢，快分享给你的朋友们吧—学习

原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让咱们一块儿并肩做战！

本公众号深耕耘生信领域多年，具备丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有须要的老师和同窗前来咨询。