對(duì)基因組測(cè)序完成后检盼，我們經(jīng)常需要統(tǒng)計(jì)測(cè)序深度（depth）和對(duì)基因組的覆蓋率（coverage）
這兩個(gè)概念有時(shí)候不太好區(qū)分氯夷，有時(shí)coverage也表示測(cè)序深度x.
我們用samtools的depth函數(shù)并結(jié)合awk來(lái)進(jìn)行統(tǒng)計(jì)赎瑰！
samtools depth -aa sample.bam >sample.coverage.txt
參數(shù)：
-a 輸出所有位點(diǎn)，包括零深度的位點(diǎn)潮酒；
-a –a义屏，--aa 完全輸出所有位點(diǎn)，包括未使用到的參考序列舰褪；
-b FILE 計(jì)算BED文件中指定位置或區(qū)域的深度皆疹；
-f FiLE 使用在FILE中的指定bam文件；
-I INT 忽略掉小于此INT值的reads占拍；
-q INT 只計(jì)算堿基質(zhì)量值大于此值的reads略就；
-Q INT 只計(jì)算比對(duì)質(zhì)量值大于此值的reads；
-r CHR:FROM –TO 只計(jì)算指定區(qū)域的reads晃酒；
第一列為chromosome表牢，第二列為position，第三列為該位點(diǎn)測(cè)序到的reads數(shù)目贝次。
有了這個(gè)信息我們就可以進(jìn)行depth和coverage分別統(tǒng)計(jì)了崔兴。
image.png
#總堿基數(shù)
wc -l sample.coverage.txt
#覆蓋到的位點(diǎn)數(shù)
awk -F "\t" 'BEGAIN{a=0}{$3>0(a++)}END{print a}' sample.coverage.txt
#覆蓋到的位點(diǎn)平均測(cè)序深度
awk -F "\t" 'BEGAIN{a=0, b=0}{$3>0(a++;b+=$3)}END{print b/a}' sample.coverage.txt
##a用于統(tǒng)計(jì)測(cè)到的位點(diǎn)，b用于統(tǒng)計(jì)測(cè)到的位點(diǎn)深度
##某條染色體（如chr1）測(cè)到的位點(diǎn)的平均測(cè)序深度
awk -F "\t" 'BEGAIN{a=0, b=0}{$3>0 && $1~/chr1$/ (a++;b+=$3)}END{print b/a}' sample.coverage.txt
能做這種分析的各種工具很多蛔翅，bedtools敲茄，samtools，GATK等山析，但是我實(shí)在記不住這么多命令折汞，就用這種簡(jiǎn)單的腳本進(jìn)行統(tǒng)計(jì)吧~~~
