【直播】我的基因组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
校对编辑:一只思考问题的熊
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 详解如何使用Spark和Scala分析Apache访问日志
- Centos 7.4中的远程访问控制的实现方法
- Linux查找特定程序whereis实例详解
- SpringBoot实战(二):SpringMvc接收xml请求
- CentOS环境下安装PHPUnit的方法分析
- Centos下安装Ansible的示例代码
- ubuntu18.04安装搜狗拼音的简易教程
- linux中$符号的基础用法总结
- Linux下的 mariadb 使用 root 用户启动方式(推荐)
- window与linux项目部署之linux文件路径不存在问题
- Ubuntu 18.04安装 pyenv、pyenv-virtualenv、virtualenv、Numpy、SciPy、Pillow、Matplotlib
- Linux使用 iftop 实时监控网卡的流量
- Centos中TCPWrappers访问控制实现
- CentOS7 重新分配分区大小的实现方法
- Linux 下安装pip包的方法