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

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

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

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

---

以下为分析过程：

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

1. 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代码如下：

1. a=read.table('../tmp.txt')

2. a$GC = a[,4]/a[,3]

3. a$depth = a[,5]/a[,3]

4. a = a[a$depth<100,]

5. plot(a$GC,a$depth)

6. library(ggplot2)

7. # GET EQUATION AND R-SQUARED AS STRING

8. # SOURCE: http://goo.gl/K4yh

9. lm_eqn <- function(x,y){

10. m <- lm(y ~ x);

11. eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,

12. list(a = format(coef(m)[1], digits = 2),

13. b = format(coef(m)[2], digits = 2),

14. r2 = format(summary(m)$r.squared, digits = 3)))

15. as.character(as.expression(eq));

16. }

17. p=ggplot(a,aes(GC,depth)) + geom_point() +

18. geom_smooth(method='lm',formula=y~x)+

19. geom_text(x = 0.5, y = 100, label = lm_eqn(a$GC , a$depth), parse = TRUE)

20. p=p+theme_set(theme_set(theme_bw(base_size=20)))

21. p=p+theme(text=element_text(face='bold'),

22. axis.text.x=element_text(angle=30,hjust=1,size =15),

23. plot.title = element_text(hjust = 0.5) ,

24. panel.grid = element_blank(),

25. #panel.border = element_blank()

26. )

27. print(p)

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

文：Jimmy

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