【直播】我的基因组51:画全基因范围内的染色体reads覆盖度图

前面我们已经详细讲解过如何根据窗口来统计每条染色体的每个片段的GC含量，还有平均测序深度，请大家自行前往前面查看脚本及实现方式！【直播】我的基因组47:测序深度和GC含量的关系

那么如果得到了如下的数据：

1. > head(dat)

2. chr number length GC counts depth

3. 1 chrY 215 98427 663 1443853 14.66928

4. 2 chr3 445 99517 17945 3906339 39.25298

5. 3 chr6 130 99698 24282 3342284 33.52408

6. 4 chr5 328 98445 16271 3071504 31.20020

7. 5 chr4 1035 99867 23631 3067318 30.71403

8. 6 chr16 582 99910 19495 3585809 35.89039

很明显，上面是以100Kb区域为窗口，我们就需要进行可视化，如下：

（抱歉，画的还是有点丑，可视化的确不是我擅长的！）

这个图有很多需要改进的地方，比如X坐标轴应该对每一个染色体来说都不一样，染色体的长度很明显可以看出来的， 但是我简单粗暴的取了最长染色体的长度！配色也不好看，大家可以参照Y叔的<食色性也>去寻找最佳配色，我实在是太忙了，没空做这些美化了。

从上面的图，我们可以得到很多信息：

1号染色体中间的测序深度有点不稳定；

9号染色体中间有一大块测序深度明显偏低，需要后面详细探究；

13,14,15,21,22号染色体开头处有大片段覆盖度为0的情况，也行是端粒处没有测到，或者这些地方就是N碱基。也需要仔细探究。

R代码如下：

1. rm(list = ls())

2. file <- 'filter-bam/GC_stat.100k.txt'

3. dat <- read.table(file, sep = "\t", fill=TRUE,stringsAsFactors = F)

4. colnames(dat)=c('chr','number','length','GC','counts')

5. keep_chr <- dat$chr %in% c(paste0('chr',c(1:22,'X','Y')))

6. dat <- dat[keep_chr,]

7. dat$depth <- dat$counts/dat$length

8. library(ggplot2)

9. head(dat)

10. # To change plot order of facet wrap,

11. # change the order of varible levels with factor()

12. dat$chr <- factor(dat$chr , levels =c(paste0('chr',c(1:22,'X','Y'))) )

13. png('coverage.png',height = 1000,width = 1000)

14. p <- ggplot(dat,aes(number,depth))+geom_area(aes(fill=chr))+ylim(0, 100)

15. p <- p+facet_wrap( ~ chr,ncol=2)

16. print(p)

17. dev.off()

上面得到的是以100Kb为窗口的统计结果，如果我们以10Kb来统计绘图，结果如下：

肉眼上，几乎看不出什么区别，同样的代码，我就不重复show啦。

（虽然我还统计了以1Kb为窗口结果，但是不想画图了，感觉都差不多了，而且1Kb的窗口统计结果文件有77Mb，画图挺耗费时间的。）

文：Jimmy

图文编辑：吃瓜群众