samtools depth 用於外顯子未覆蓋區域的統計及統計未覆蓋區域的意義
samtools depth主要用來從bam檔案中統計指定區域的深度情況。
首先還是簡單介紹一下samtools depth的基本用法,如下圖所示
我們可以通過samtools depth option 1.bam 2.bam...的方式來執行該軟體,此外,最常用的引數是-r引數,我們可以指定一些區域來生成指定區域的深度情況,也可以通過輸入一個-b引數輸入一個bed 檔案來實現該過程。
而得到的結果一共由3列構成,chrome,postion和深度,如下圖所示:
實際工作中遇到的一個問題是使用samtools depth來獲取捕獲區域的未覆蓋區域,思路如下:
- 使用samtools depth -a 引數來獲取捕獲bed區域的缺失區域
- 使用samtools depth -aa 引數來獲取捕獲bed區域的覆蓋度為0的區域(包括未覆蓋和缺失區域)
- 使用覆蓋度為0的區域減去缺失區域得到未覆蓋區域
另外一個問題是捕獲效率值的由來,在測序實驗開始的時候,會使用捕獲試劑盒來捕獲想要的區域的片段,但是這些片段並不是就肯定在想要的區域內的,由於序列的相似性也會被捕獲到,最終我們測序得到的reads,將這些reads比對回參考基因組就發現有很多reads是比對不到捕獲區域的,而使用比對得到的所有鹼基數比上比對到捕獲區域的鹼基數就是捕獲效率值。
這裡我們可以來舉一個提取未覆蓋區域的例子:
我們找到一個未覆蓋區域,然後分別使用samtools depth -a和-aa引數去試一下
samtools-1.6 depth -a -r 20:2673783-2673784 test.final.bam
#無結果
samtools-1.6 depth -aa -r 20:2673783-2673784 test.final.bam
#結果
20 2673783 0
20 2673784 0
再找一個純合缺失的位置去看一下:
samtools-1.6 depth -aa -r 1:977201-977202 test.final.bam 1 977201 0 1 977202 0 samtools-1.6 depth -a -r 1:977201-977202 test.final.bam 1 977201 0 1 977202 0
發現針對純合缺失位點,兩種方式都可以檢出該區域的深度。
最後一個問題是為什麼要統計未覆蓋區域呢?這是因為全外顯子測序對於一些大片段的缺失的檢測能力是有限的,有時候,一整個基因的缺失,或者大於1kb的缺失往往會漏掉,通過統計uncover區域的情況,我們可以檢測是否有這種大片段的缺失,通過使用phenolyzer來對uncover區域的基因註釋進行關聯,很可能發現這一類的大片段缺失,有實際的例子,就是通過這種方式發現了一個重要的關聯致病基因的缺失