测序深度分布图及代码（shell和R）

微信公众号：生物信息学习

有一个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

测序深度的分布图就出来了。

微信公众号：生物信息学习