【直播】我的基因组(20):覆盖度详细探究

前面我们在第8讲提到了公司给我的一个报告的统计表格,有人反映不会做。

本来应该只需要给我6亿条reads的(PE150测序,人30X),但是足足给了我8.9亿条!(但事实上很多paper发表的基因组高于60X的也不少)

表格里面提到了好几个概念,比如duplicate的reads,一般说的PCR造成的duplicate,在找变异的时候需要去除掉。然后是那些比对到了不同染色体的reads pair,虽然只有2.29% ,也是需要重点分析的。(前面我也讲了如何提取以及分别分析它们!)

如果只是想重新前面的这些统计指标,非常很简单,就是samtools工具提供了一个flagstat功能,用法及结果如下:

samtools flagstat P_jmzeng.final.bam

 

899361748 + 0 in total (QC-passed reads + QC-failed reads)

8597742 + 0 secondary

0 + 0 supplementary

132556557 + 0 duplicates

890858540 + 0 mapped (99.05% : N/A)

890764006 + 0 paired in sequencing

445382003 + 0 read1

445382003 + 0 read2

853255862 + 0 properly paired (95.79% : N/A)

881249604 + 0 with itself and mate mapped

1011194 + 0 singletons (0.11% : N/A)

20382234 + 0 with mate mapped to a different chr

12511988 + 0 with mate mapped to a different chr (mapQ>=5)

从结果很明显可以看到,公司的确给了我8.9亿的reads,这个没错,duplicate的情况也的确是132556557/899361748 = 14.73%,而mapping的情况和proper mapping的情况也都显示好了,可以看到公司用的是同样的命令和方法!没有什么神秘的,我们生信工程师,做得就是这个,而且可以做得更好。

后面的探索全基因组区域中碱基覆盖深度不低于多少X的比例,是需要画一个图,有非常多的现成的工具可以使用,包括 BedTools' genomeCov 、 GATK's DepthOfCoverage,还有Picard suite的几个命令。

其实本身原理很简单,就是把全基因组的每个坐标的depth都得到,然后得到depth的频数,然后画图。

我们可以对每条染色体单独来绘图,也可以针对全基因组来绘图,当然,公司给我们的数据是针对全基因组的。 脚本很简单:

samtools mpileup P_jmzeng.final.bam |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' 当是需要运行很久,毕竟全基因组的bam文件太大了。

如果分开运行,可以对下面的各条染色体bam文件批量跑一个脚本

脚本如下:

ls P_jmzeng.final.REF*.bam |while read id

do

echo $id

samtools mpileup $id |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' >$id.depth.txt

done

跑完之后,对每条染色体都会输出如下文件:

over100 110789

0 27065

1 1730286

2 2219409

3 2728526

4 3251046

5 3774335

6 4303971

~~~~~~~~~~~~~~~~~~~~~~~

后面省略94行,每一行都是两列,第一列是测序深度,第二列是有着该测序深度的位点是多少个!所以行的第二列加起来,就是染色体的长度!!!面这个例子是X染色体的,绘图如下,可以看到X染色体的测序深度其实并不怎么好,全基因组测序深度平均高达44X,可是这个X染色体超过44X的只占极少的比例。

我用excel表格简单画个图如下:(当然,作为一个高级生物信息学工程师,用excel表格是有点low,但是这里只是为了说明一个问题,我们后面还是写程序的,用R语言)

比如我们再看看10号染色体,我随意在R里面画了个图:

a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)

plot(a[,2],type = 'l',xaxt="n",lwd=3)

axis(1, at = 1:nrow(a),labels =a[,1] , las=1)

abline(v=44)

可以看到10号染色体的测序深度要显著好于X染色体,大部分的测序深度都超过了30X!!(不过,这里虽然用来R,但是出图很丑,虽然不是专业做可视化的,但是想调整美观一点问题也不大,就是好耗费时间。)

做出depth的累积曲线图了,也是很简单的!

a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)

over100=a[1,]

a=a[-1,]

a=rbind(a,over100)

a[,3]= cumsum(a[,2])

a[,4]= 1-a[,3]/sum(a[,2])

plot(a[,4],type = 'l',xaxt="n",lwd=3)

axis(1, at = 1:nrow(a),labels =a[,1] , las=1)

abline(v=10);abline(v=20);abline(v=30)

最后得到的数据如下:第4列,就是不低于多少X的比例

可以看到大于10X测序深度的比例仍然高达92.88%,效果杠杠的!!


这个测序深度累积分布图,就是很多人的重点!!!!(请大家仔细学习我上面的统计脚本和画图方法)

加油吧,骚年!!!!

请扫描以下二维码关注我们,获取直播系列的所有帖子!

菜鸟团公众号二维码

菜鸟团公众号二维码

参考:

https://www.biostars.org/p/104063/

http://www.gettinggeneticsdone.com/2014/03/visualize-coverage-exome-targeted-ngs-bedtools.html

 

Comments are closed.