【直播】我的基因组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

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-24

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Young Dreamer

针对iPhone的pt、Android的dp、HTML的css像素与dpr、设计尺寸和物理像素的浅分析

  最近被一朋友问到:css中设置一DOM的height:65px,请问显示的高度是否和Android的65dp的元素等高?脑子里瞬间闪现了一堆的概念,如dpr...

1955
来自专栏Android群英传

贝塞尔Loading——化学风暴

451
来自专栏欧阳大哥的轮子

地图开发知识之-投影坐标

由于地球是一个赤道略宽两极略扁的不规则的梨形球体,表面是一个不可展平的曲面,而地图通常是二维平面,因此在地图制图时首先要考虑把曲面转化成平面。然而,从几何意义上...

713
来自专栏ThoughtWorks

前端不止:Retina屏幕下两倍图

所见不一定即所得 眼睛是心灵的窗户,也是蒙蔽你的一种途径。 假设,我给你一张图片,你觉得肉眼可以观察到全部的细节吗? ? ---- 屏幕上一张清晰的图片 肉眼...

3455
来自专栏大数据钻研

canvas图形绘制之星空、噪点与烟雾效果

一、三合一 三个效果合成一篇文章。 有多个小伙伴问我,为何不开个公众号,现在都是移动时代,你博客文章写好后,公众号再复制一份,花不了多长时间,同时传播方便迅速,...

2944
来自专栏诸葛青云的专栏

迪杰斯特拉(dijkstra)c语言实现方法

迪杰斯特拉(dijkstra)是用来实现查找一个点到其它点最短路径的一种方法。通过查找从起点到最短距离的点,然后将该点放入到集合中,代表以及找到起点到这一点的最...

612
来自专栏MixLab科技+设计实验室

设计师会编程、程序员懂艺术:Semi Flat Design

本公众号定期更新关于 设计师、程序员发挥创意 互相融合的指南、作品。 主要技术栈: nodejs、react native、electron 本号正在更新的系列...

2946
来自专栏数据小魔方

R预设配色系统及自定义色板

关于配色的话题,已经聊过很多次了,但是就像是之前说过的,对于图形可视化而言,配色决定着作品的“颜值”,谈再多次都不嫌多。 今天是R语言配色系统综合篇的上篇(当然...

4379
来自专栏算法+

双边滤波算法的简易实现bilateralFilter

没怎么看过双边滤波的具体思路,动手写一写,看看能不能突破一下。 最后,感觉算法还是要分开 水平 与 垂直 方向进行分别处理,才能把速度提上去。 没耐性写下去了,...

3996
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 1

基因组浏览器是高通量测序分析的一个重要的可视化工具。相比于最终提供的表格,基因组浏览器可以提供更多的信息,如直观展示突变位点、查看有无新转录本或新的可变剪接形式...

832

扫码关注云+社区