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

相关文章

来自专栏月色的自留地

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

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

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

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

3646
来自专栏Y大宽

cBioportal中文教程

大规模的癌症基因组计划,比如The cancer genome atlas(TCGA) and the International cancer genome ...

2552
来自专栏企鹅号快讯

麦子陪你做作业(二):KEGG通路数据库的正确打开姿势

KEGG是通路数据库中最庞大的,涵盖基因组网络信息,主要注释基因的功能和调控关系。当我们选到了合适的候选分子,单变量研究也已做完,接着研究机制的时便可使用到它。...

2726
来自专栏机器之心

资源 | Facebook开源DrQA的PyTorch实现:基于维基百科的问答系统

选自GitHub 机器之心编译 参与:Panda 今年 4 月,斯坦福大学和 Facebook 人工智能研究所在 arXiv 发布了一个基于维基百科的开放域问题...

3785
来自专栏企鹅号快讯

麦子陪你做作业(二):KEGG通路数据库的正确打开姿势

KEGG是通路数据库中最庞大的,涵盖基因组网络信息,主要注释基因的功能和调控关系。当我们选到了合适的候选分子,单变量研究也已做完,接着研究机制的时便可使用到它。...

2670
来自专栏数据小魔方

R语言可视化——数据地图应用(东三省)

今天是一个案例应用,采用东北三省地图进行离散颜色映射,让大家感受下R语言在地理信息空间可视化方面的强大功能,同时也会对之前强调过的地图配色技巧进行应用。 加载工...

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

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

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

4408
来自专栏生信技能树

单细胞转录组3大R包之monocle2

主要是针对单细胞转录组测序数据开发的,用来找不同细胞类型或者不同细胞状态的差异表达基因。分析起始是表达矩阵,作者推荐用比较老旧的Tophat+Cufflinks...

2.5K11
来自专栏龙行天下CSIEM

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

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

1522

扫码关注云+社区