如何从BAM文件绘制ATAC-seq插入片段长度分布图？

欢迎关注“生物信息学实践指南”！在ATAC-seq（Assay for Transposase-Accessible Chromatin using sequencing）数据处理过程中，插入片段长度的分布图能够帮助研究者评估实验的质量。这种图表通常显示出每200碱基对(bp)出现一个峰值，这一周期性的模式反映了核小体的数量分布。

在进行数据分析时，首先需要从BAM文件中提取插入片段的长度信息。对于双端测序数据，BAM文件的第9列包含了插入片段的长度。提取这些信息后，可以通过R语言等工具绘制出插入片段长度的分布图。

下面是一段用于提取插入片段长度的命令：

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

这段命令首先使用samtools从BAM文件中读取比对信息，然后通过awk脚本计算每个读段的绝对插入片段长度，并去重以避免同一片段的多次计数。最终结果保存在fragment.length.txt文件中。

接下来，使用R语言绘制插入片段长度的分布图：

data <- read.table("fragment.length.txt", header = FALSE)
# 设置插入片段长度的阈值，过滤掉过长的片段
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 Distribution")

此外，还可以使用Picard的CollectInsertSizeMetrics或Bedtools的bamPEFragmentSize工具来计算插入片段长度，但上述方法更为直接且易于实现。

更多生物信息学相关的内容，请持续关注我们的公众号，获取最新资讯和技术分享。