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

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏量子位

秘籍:如何用廉价硬件玩转深度学习,成本不到1000美元

作者Lukas Biewald,是CrowdFlower创始人。 量子位编译整理。 问:搭建一个深度学习系统拢共要花多少钱? 答:在树莓派上运行TensorFl...

19810
来自专栏生信技能树

一篇文章学会miRNA-seq分析

第一讲:文献选择与解读 前阵子逛BioStar论坛的时候看到了一个关于miRNA分析的问题,提问者从NCBI的SRA中下载文献提供的原始数据,然后处理的时候出现...

6626
来自专栏AI研习社

用机器学习搞艺术,谷歌 Megenta 项目集锦(附 Github)

雷锋网 AI 研习社按:本文为雷锋字幕组编译的技术博客,原文 Make Music and Art Using Machine Learning,作者 mage...

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

近期GitHub上最热门的开源项目(附链接)

来源:开源最前线 2 月份 GitHub 上最热门的开源项目又出炉了,又有哪些新的项目挤进热门榜单了呢,一起来看看。 ……………………………… 1、nocode...

3509
来自专栏安全领域

在物联网中应用机器学习:使用 Android Things 与 TensorFlow

在本教程中,我们将探索如何使用 Android Things 和 TensorFlow 将机器学习应用到物联网中。

71816
来自专栏小詹同学

20秒画完小猪佩奇“社会人”,程序猿的手法是你想不到的独特

今年社交平台上最火的带货女王是谁?范冰冰?杨幂?Angelababy?不,是猪猪女孩小猪佩奇。

1081
来自专栏Java架构师学习

一些最好用的数据可视化工具

摘要: 如今同质化的应用越来越多,应用开发者也开始在用户体验上下功夫,比如数据可视化,将一大堆密密麻麻的数字转成图表形式,可以更直观地向用户展示数据之间的联系和...

4985
来自专栏IT派

7月Python最佳开源项目Top 10

【导读】七月就要结束了,小编为大家整理了本月 Python 最受欢迎的十大开源项目。他山之石,可以攻玉,爱好Python的朋友们一起学习Github上的优秀项目...

573
来自专栏JackeyGao的博客

2016年Python十大文章

在过去一年, 我们对10000篇Python相关的文章进行了排名, 并选择出排名前十的文章. (0.1%的几率), 可以帮助您在2017年推进你的技术生涯.

511
来自专栏木子昭的博客

泰坦尼克乘客存活状况(决策树案例)

1912年4月15日凌晨2点20分,“永不沉没”的“泰坦尼克”走完了它短暂的航程,缓缓沉入大西洋这座安静冰冷的坟墓。 ? 欢迎你们说我幼稚荒诞,也欢迎你...

34312

扫码关注云+社区