对基因组测序完成后，我们经常需要统计测序深度（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等，但是我实在记不住这么多命令，就用这种简单的脚本进行统计吧~~~