搜索
查看: 4329|回复: 1

【直播】我的基因组19:根据比对结果来统计测序深度和覆...

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-2 13:20:40 | 显示全部楼层 |阅读模式
【直播】我的基因组19:根据比对结果来统计测序深度和覆盖度
看来本次直播我的基因组分析流程效果还不错,不少朋友跟着自己动手开始分析全基因组测试数据了,值得表扬。其中有好几个朋友留言向我反映公司给的统计报告里面的覆盖度的问题,如下图:


我在第 8讲中写道,每条染色体的覆盖度都接近于100%,而且测序深度大多在40X以上,不少读者表示看晕了,明明这个条形图显示覆盖度不到60%呀!其实是公司的这个图没有做好,它里面的覆盖度用的是最上面的曲线显示的,条形图是测序深度!!!
测序深度和覆盖度的示意图如下:


那么我们就来自己动手统计一下比对好的sam/bam文件的测序深度和覆盖度吧!
这个统计主要依赖于samtools的depth功能,或者说mpileup功能,输入文件都是sort好bam格式的比对文件。事实上,你如果读samtools的源代码就会发现,其实depth功能调用的就是mpileup的函数。但是mpileup可以设置一系列的过滤参数。而depth命令是纯天然的,所以mpileup的结果一定会小于depth的测序深度。对mpileup,可以不选择-u -f 参数指定参考基因组,因为我们只需要测序深度情况,还有,可以指定-q 1 来过滤掉多比对情况。还可以用-Q来过滤低质量的碱基(base pair),用-A来过滤无法定位的reads,结果如下:


针对这个全基因组位点的统计结果,我们很容易写脚本来计算每条平均的测序深度和各个染色体的覆盖度。

[Perl] 纯文本查看 复制代码
nohup time samtools mpileup P_jmzeng.final.bam |perl -alne '{$pos{$F[0]}++;$depth{$F[0]}+=$F[3]} END{print "$_\t$pos{$_}\t$depth{$_}" foreach sort keys %pos}' 1>coverage_depth.txt 2>coverage_depth.log &

nohup ... & ...为命令 表明命令后台执行  也可以
nohup ... & > *.log 将运行结果 写入到一个log文件里面

time 命令 可以统计 命令 运行的时间

这个脚本运行会比较慢,因为是针对整个55G的bam文件。耗时如下:
real        130m34.855s
user        237m14.308s
sys        1m45.692s

结果如下:


其中chromosome的长度在bam文件里面可以看到,用samtools view -H P_jmzeng.final.bam 即可!!!把上面的表格可视化,就是文章最开头的figure。但是很明显公司给我的各个指标均高于我自己算出来的!尤其是其中几个染色体的coverage非常差。当然平均测序深度,我在这里选择的是整个染色体的长度作为分母,对于那些覆盖度很差的染色体,平均测序深度就被被拉低。所以我起初猜想是不是问题出在我用的是samtools的mpileup命令,而不是depth命令!(因为我以前从来没有如此细致的去比较它们的差别,其实这这个命令的确有差别,但是对全基因组层面的统计指标几乎没什么影响)

下面的表格,是我用depth命令的结果,而且平均测序深度我选择被覆盖到的染色体长度作为分母,所以看到测序深度有些许提升,但是有几条染色体的覆盖度堪忧。与公司给我的报告不符合!


先不要急着怪公司,请听下回分解

这个问题我会详细讲解,请关注后面的帖子:
或者不用重复造轮子,bedtools有这个命令:[url=]http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html[/url]
而且Qualimap软件也可以做到,后面我们会讲到!

参考:
[url=]http://seqanswers.com/forums/showthread.php?t=17438[/url]
[url=]https://www.biostars.org/p/218049/[/url]
[url=]https://www.biostars.org/p/3326/[/url]
[url=]http://www.danielecook.com/calculate-depth-coverage-bam-file/[/url]
[url=]Should Samtools Pileup Be Performed On Uniquely Mapped Reads Or All The Reads?[/url]
[url=]Genomic Coverage - Samtool's undocumented "depth" verses the poorly documented pileup.[/url]
[url=]Discrepancy In Samtools Mpileup/Depth And Bedtools Genomecoveragebed Counts[/url]


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:R语言初学笔记之吾日三省吾身!
下一篇:【直播】我的基因组20:覆盖度详细探究
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

13

主题

44

帖子

522

积分

高级会员

Rank: 4

积分
522
发表于 2018-4-16 17:54:05 | 显示全部楼层
为什么我用这个脚本跑出来只有三列啊,不知道这三列是什么意思
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-20 19:27 , Processed in 0.030465 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.