【直播】我的基因组47:测序深度和GC含量的关系

在前面我们用 ChIP-seq 的分析方法可视化了一下我的 WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关

比如在TSS附近,就很明显看到了一个测序深度峰值(具体内容点击 【直播】我的基因组 44:比对文件画profile和heatmap图),但是前面我们并没有给出直接的解答而是简单的提到这是二代测序的特点——GC含量片段偏好性

作为一个合格的生物信息学工程师,我当然要把这个理论用自己的代码和数据来亲身实践一遍。

以下为分析过程:

首先,把全基因组的bam文件用 mpileup模式输出,根据 1000bp 的窗口滑动来统计每个窗口的测到的碱基数,GC碱基数,测序总深度!(代码比较复杂,一般人可能理解不来)

samtools mpileup -f ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta ../P_jmzeng.final.bam|head -1000000 |perl -alne '{$pos=int($F[1]/1000); $key="$F[0]\t$pos";$GC{$key}++ if $F[2]=~/[GC]/;$counts_sum{$key}+=$F[3];$number{$key}++;}END{print "$_\t$number{$_}\t$GC{$_}\t$counts_sum{$_}" foreach sort{$a<=>$b} keys %number}'

上面的代码写的不好,跑10万行需要 4s,跑一百万行需要36s,我估计把这8.9亿行的bam运行完,这样推算是10小时即可,但事实上我已经跑了一整天了!我感觉自己的脚本能力在面对大数据(300Gb的全基因组)有点捉鸡!

不过不要紧,我们就拿前面的百万行数据做一个测试就好了。

结果如下:

说明 前面两行是窗口的坐标,第几号染色体的第几个窗口,后面3行是数据,分别是每个窗口的测到的碱基数,GC碱基数,测序总深度。

接下来,将上面的文件导入到 R里进行可视化。

PS:这个线性回归图不会看的,自己去搜索或者去看生信技能树论坛的文章: http://www.biotrainee.com/thread-695-1-1.html (复制链接到浏览器打开或者点击最下方的阅读原文)。 我觉得我这次画的图还不错,很明显能看到这个趋势,GC含量比较高的窗口,有着相应比较高的测序深度!

至此,完美的证明了文章开头的结论!

给自己一百个赞,虽然我没有对全基因组数据做验证,但是基因组差异并没有很大,我也随机抽样测试了几次都有这个趋势。

最后,给出我的 R代码如下:

a=read.table('../tmp.txt')a$GC = a[,4]/a[,3]a$depth = a[,5]/a[,3]a = a[a$depth<100,]plot(a$GC,a$depth)library(ggplot2)# GET EQUATION AND R-SQUARED AS STRING# SOURCE: http://goo.gl/K4yhlm_eqn <- function(x,y){m <- lm(y ~ x);eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,list(a = format(coef(m)[1], digits = 2),b = format(coef(m)[2], digits = 2),r2 = format(summary(m)$r.squared, digits = 3)))as.character(as.expression(eq));}p=ggplot(a,aes(GC,depth)) + geom_point() +geom_smooth(method='lm',formula=y~x)+geom_text(x = 0.5, y = 100, label = lm_eqn(a$GC , a$depth), parse = TRUE)p=p+theme_set(theme_set(theme_bw(base_size=20)))p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1,size =15),plot.title = element_text(hjust = 0.5) ,panel.grid = element_blank(),#panel.border = element_blank())print(p)

关于画图,大家可以参考下面这个链接:http://stackoverflow.com/questions/7549694/ggplot2-adding-regression-line-equation-and-r2-on-graph

文:Jimmy

校对编辑:一只思考问题的熊

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

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

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信宝典

生信宝典之傻瓜式 (二) 如何快速查找指定基因的调控网络

我是谁?我在哪儿?我在查什么? 在信息爆炸的时代,相信很多小伙伴在查文章时会因信息量太大而抓狂。今天带来一款设计简洁、功能全面的基因功能查询工具,助你事半功倍,...

1866
来自专栏phodal

如何为技术博客设计一个推荐系统(中):基于 Google 搜索的半自动推荐

与统计学相比,基于内容来向用户推荐相似的内容,往往更容易获得。对于推荐来说,则有两种方式: 手动推荐 自动推荐 (PS:我承认,这句话说了等于没说。) 如下图所...

2016
来自专栏目标检测和深度学习

学者必备5大科研神器

觉得不错的话 欢迎转发 点赞 收藏 留言 ? 01 MathType ? 一款强大的数学公式编辑器 推荐指数:⭐⭐⭐⭐⭐ 写论文的时候 最烦的就是插入公式! ...

3268
来自专栏钱塘大数据

钱塘干货 | 数据收集和处理工具一览

进入大数据时代,调查报道愈加成为信息战。从哪里收集有效数据?如何抽取、筛选、整合、分类大量琐碎的信息?如何分享、存储数据,并实现随取随用?钱塘君整理了一张数据收...

3677
来自专栏分子生物和分子模拟计算

荧光探针分子与蛋白质的作用

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

【学习】Python可视化工具概述-外文编译

本文由 PPV课 - korobas 翻译,未经许可,禁止转载! 原文翻译链接:http://pbpython.com/visualization-tools...

3687
来自专栏Data Analysis & Viz

Gephi绘制微博转发图谱:以“@老婆孩子在天堂”为例

以前看过一篇提取《釜山行》剧本中人物,并用Gephi绘制关系图谱的文章,因此想用Gephi绘制下微博转发情况,借此来换个角度看看微博内容是怎么扩散的。其中爬取转...

1052

Python机器学习的生态系统

Python生态系统正在不断成长,并可能成为机器学习的统治平台。

1867
来自专栏理论坞

尼尔森十大可用性原则知多少?

尼尔森(Jakob Nielsen)是一位人机交互学博士(Technical University of Denmark in Copenhagen), 于19...

773
来自专栏数据科学与人工智能

【陆勤践行】Python和数据科学的起步指南

Python拥有着极其丰富且稳定的数据科学工具环境。遗憾的是,对不了解的人来说这个环境犹如丛林一般(cue snake joke)。在这篇文章中,我会一步一步指...

21310

扫描关注云+社区