对基因组测序完成后，我们经常需要统计测序深度（depth）和对基因组的覆盖率（coverage）
这两个概念有时候不太好区分，有时coverage也表示测序深度x.
我们用samtools的depth函数并结合awk来进行统计！
samtools depth -aa sample.bam >sample.coverage.txt
参数：
-a 输出所有位点，包括零深度的位点；
-a –a，--aa 完全输出所有位点，包括未使用到的参考序列；
-b FILE 计算BED文件中指定位置或区域的深度；
-f FiLE 使用在FILE中的指定bam文件；
-I INT 忽略掉小于此INT值的reads；
-q INT 只计算碱基质量值大于此值的reads；
-Q INT 只计算比对质量值大于此值的reads；
-r CHR:FROM –TO 只计算指定区域的reads；
第一列为chromosome，第二列为position，第三列为该位点测序到的reads数目。
有了这个信息我们就可以进行depth和coverage分别统计了。
image.png
#总碱基数
wc -l sample.coverage.txt
#覆盖到的位点数
awk -F "\t" 'BEGAIN{a=0}{$3>0(a++)}END{print a}' sample.coverage.txt
#覆盖到的位点平均测序深度
awk -F "\t" 'BEGAIN{a=0, b=0}{$3>0(a++;b+=$3)}END{print b/a}' sample.coverage.txt
##a用于统计测到的位点，b用于统计测到的位点深度
##某条染色体（如chr1）测到的位点的平均测序深度
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等，但是我实在记不住这么多命令，就用这种简单的脚本进行统计吧~~~
