微信公众号：生物信息学习
有一个bed文件，想要知道这个bed文件对应的测序深度的分布，bed文件如下图所示：
chr1 4427277 4427583
chr1 4429216 4429458
用samtools对这个区间内的每个点做深度计算，代码如下：
samtools depth test_sorted_dedup.bam -b test.bed >test.stat
结果如下所示（其中第一列为染色体，第二列为position，第三列即为这个position的深度）：
chr10 3787066 109
chr10 3787067 108
chr10 3787068 109
chr10 3787069 111
chr10 3787070 111
chr10 3787071 111
chr10 3787072 111
chr10 3787073 112
chr10 3787074 114
chr10 3787075 115
chr10 3787076 116
chr10 3787077 118
用shell对区间内所有的测序深度的depth做统计，代码如下：
cut -f3 test.stat |sort|uniq -c >test.stat.count
统计后的结果如下所示：
3503 55
5097 1
2508 265
3607 42
4453 8
4694 7
1300 2762
2898 165
其中第二列为depth，第一列为这个depth的频数，如第一行，即表示有3503个点的深度为55.
接下来用R做这些深度分布图，代码如下：
library(ggplot2)
mytest<-read.table("test.stat.count")
names(mytest)<-c("depth","count")
ggplot(mytest,aes(x=depth,y=count))+geom_point(col="red")
结果如下所示：
image
测序深度的分布图就出来了。
