|
計算測序深度是數(shù)據(jù)分析中的一個常規(guī)操作，常用的工具有以下幾種 1. samtools depth命令如下 samtools depth input.bam > sample.depth2. samtools mpileup命令如下 samtools mpileup -A -Q 0 chrY.bam | cut -f1,2,4 > sample.depth3. bdetools genomecov命令如下 bedtools genomecov -ibam input.bam -d > sample.depth盡管方法不同,，但是輸出結(jié)果的格式是相同的,，示意如下 第一列是染色體名稱，第二列是染色體上的坐標(biāo),，第三列是對應(yīng)的測序深度,。原本以為計算測序深度就是這么輕松的一件事，但是在比較不同方法的輸出結(jié)果時,，卻發(fā)現(xiàn)部分區(qū)域samtools計算的結(jié)果和bedtools的結(jié)果對應(yīng)不上,，結(jié)果如下 第三列為samtools軟件計算出來的結(jié)果，第四列為bedtools軟件計算出來的結(jié)果,，可以看到,，samtools的結(jié)果比bedtools少了很多,。 為了確定這段區(qū)域真實的測序深度，將bam文件導(dǎo)入IGV軟件之中進行查看,，對應(yīng)區(qū)域的結(jié)果如下 從reads來看,，確實應(yīng)該是10和16條，那么samtools計算出來的結(jié)果又是什么,， 總不可能是samtools的bug吧,，作為一個應(yīng)用這個廣泛的軟件，不可能有如此低級的錯誤,。 從結(jié)果來看,，samtools在計算depth的過程中對reads進行了過濾，那么它過濾的原則是什么呢,？通過查看bam文件中第二列的flag我找到了規(guī)律,，143-150bp的reads有以下10條 第二列的flag含義如下 0x91 145 PAIRED,REVERSE,READ2可以看到，其中有一個SECONDARY比對,，對應(yīng)數(shù)字355,，這樣的reads有4條，去除這4條之后,，就是samtools計算的最終結(jié)果了,。 作為SAM文件的官方操作工具，samtools內(nèi)置了對flag的過濾,， 而bedtools默認沒有進行這樣的過濾,，只是簡單統(tǒng)計了該區(qū)域比對上的reads數(shù)目。 在使用IGV展示測序深度時,，可以用igvtools先計算出tdf文件再導(dǎo)入IGV中,，這樣比直接導(dǎo)入bam文件快很多。在tdf文件中,，其測序深度和bedtools軟件的計算結(jié)果是一致的,，也是只統(tǒng)計了read數(shù)目，沒有做過濾,。 ·end· |
|
|
