歡迎關注」生信修煉手冊」!
對於ATAC文庫而言，其插入片段的長度分布有著非常典型的規律，示意如下
每200bp會存在一個峰，這個周期性波動反應的是核小體的個數。在ATAC_seq的數據分析中，會對插入片段長度分布進行可視化，觀察其是否符合這樣的周期性規律，一定程度可以反映文庫構建的質量，那麼如何在做這樣一張分布圖呢?
比對之後我們會得到bam文件，畫圖所需的插入片段長度就需要從bam文件中提取，需要注意，這裡的插入片段是文庫中adapter之間的插入片段，即fragment， 需要和insert size區別開來。
對於單端測序而言，只有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代碼如下
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")
輸出結果示意如下
這種是最簡單的方式，除此之外，還有picard的CollectInsertSizeMetrics， bedtools的bamPEFragmentSize也都可以計算插入片段長度，但是在我看來，還是本文介紹的方式最為簡單直接。
原創不易，歡迎收藏，點讚，轉發！生信知識浩瀚如海，在生信學習的道路上，讓我們一起並肩作戰！本公眾號深耕耘生信領域多年，具有豐富的數據分析經驗，致力於提供真正有價值的數據分析服務，擅長個性化分析，歡迎有需要的老師和同學前來諮詢。轉發本文至朋友圈，後臺私信截圖即可加入生信交流群，和小夥伴一起學習交流。
掃描下方二維碼，關注我們，解鎖更多精彩內容！