搜索
查看: 3654|回复: 2

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-25 14:32:41 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-2-25 14:40 编辑

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

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


那么如果得到了如下的数据:
  • [AppleScript] 纯文本查看 复制代码
    > head(dat)
    chr number length GC counts depth
    1 chrY 215 98427 663 1443853 14.66928
    2 chr3 445 99517 17945 3906339 39.25298
    3 chr6 130 99698 24282 3342284 33.52408
    4 chr5 328 98445 16271 3071504 31.20020
    5 chr4 1035 99867 23631 3067318 30.71403
    6 chr16 582 99910 19495 3585809 35.89039


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


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


从上面的图,我们可以得到很多信息:
1号染色体中间的测序深度有点不稳定;

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

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

R代码如下:
  • [AppleScript] 纯文本查看 复制代码
    rm(list = ls())
    file <- 'filter-bam/GC_stat.100k.txt'
    dat <- read.table(file, sep = "\t", fill=TRUE,stringsAsFactors = F)
    colnames(dat)=c('chr','number','length','GC','counts')
    keep_chr <- dat$chr %in% c(paste0('chr',c(1:22,'X','Y')))
    dat <- dat[keep_chr,]
    dat$depth <- dat$counts/dat$length
    library(ggplot2)
    head(dat)
    # To change plot order of facet wrap,
    # change the order of varible levels with factor()
    dat$chr <- factor(dat$chr , levels =c(paste0('chr',c(1:22,'X','Y'))) )
    png('coverage.png',height = 1000,width = 1000)
    p <- ggplot(dat,aes(number,depth))+geom_area(aes(fill=chr))+ylim(0, 100)
    p <- p+facet_wrap( ~ chr,ncol=2)
    print(p)
    dev.off()


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




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

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



文:Jimmy

图文编辑:吃瓜群众



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:三代测序建库-2 : Template preparation workflow
下一篇:【直播】我的基因组52:X和Y染色体的同源区域探索
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

3

帖子

75

积分

注册会员

Rank: 2

积分
75
发表于 2017-3-16 21:26:48 | 显示全部楼层
回复

使用道具 举报

7

主题

26

帖子

759

积分

版主

Rank: 7Rank: 7Rank: 7

积分
759
发表于 2018-3-19 12:51:36 | 显示全部楼层
输入文件是什么,感觉有点奇怪
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 00:18 , Processed in 0.032874 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.