給你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·工具

—若是喜歡,快分享給你的朋友們吧—學習



原創不易,歡迎收藏,點贊,轉發!生信知識浩瀚如海,在生信學習的道路上,讓咱們一塊兒並肩做戰!
本公衆號深耕耘生信領域多年,具備豐富的數據分析經驗,致力於提供真正有價值的數據分析服務,擅長個性化分析,歡迎有須要的老師和同窗前來諮詢。
  更多精彩



  寫在最後

轉發本文至朋友圈,後臺私信截圖便可加入生信交流羣,和小夥伴一塊兒學習交流。


掃描下方二維碼,關注咱們,解鎖更多精彩內容!


一個只分享乾貨的

生信衆號




本文分享自微信公衆號 - 生信修煉手冊(shengxinxiulian)。
若有侵權,請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」,歡迎正在閱讀的你也加入,一塊兒分享。

相關文章
相關標籤/搜索