【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏ThoughtWorks

TW洞见〡3D打印的各种问题及解决方案

文章作者来自ThoughtWorks:贺思聪 ,图片来自网络。 3D打印机已经买回来几个月了,基本上每天都要打印一些东西,期间遇到了很多的问题积累了很多的经验...

38012
来自专栏月色的自留地

就算会用python画颗心,可你依然还是只单身狗

2388
来自专栏大数据挖掘DT机器学习

用python对汽车油耗进行数据分析

- 从http://fueleconomy.gov/geg/epadata/vehicles.csv.zip 下载汽车油耗数据集并解压 - 进入jupyt...

3726
来自专栏Y大宽

ToppGene Suite中文使用指南

2007.12:Improved human disease candidate gene prioritization using mouse phenoty...

1523
来自专栏生信技能树

【直播】我的基因组48:我可能测了一个假的全基因组

背景知识 男性只有一条X染色体和一条Y染色体,所以,理论上它们上面的SNV都应该是纯合的! X,Y除了同源区域外,其它地方差异很大。所以在女性样本里面即使是混入...

31712
来自专栏CreateAMind

ray框架及ray-rllab

1582
来自专栏龙行天下CSIEM

科学瞎想系列之八十四 永磁电机(6)

【图片部分来自网络如有侵权敬请邮箱联系。欢迎原文转发到朋友圈,未经许可的媒体平台谢绝转载,如需转载或合作请邮件联系。联系邮箱laolicsiem@126.com...

2092
来自专栏生信技能树

使用GSVA方法计算某基因集在各个样本的表现

算法本身就不是很好理解,并不强求一定要理解透彻,可以参考2005年的GSEA算法:

2853
来自专栏生信技能树

比对到hg19和hg38对somatic变异的寻找影响很大

其中B是正常组织的WES数据,使用varscan找somatic mutation的时候作为normal,然后对另外两个样本(D和T)计算。 从这个bam文件可...

1723
来自专栏PPV课数据科学社区

用python对汽车油耗进行数据分析

- 从http://fueleconomy.gov/geg/epadata/vehicles.csv.zip 下载汽车油耗数据集并解压 - 进入jupyte...

4538

扫码关注云+社区