【直播】我的基因组77：批量计算每个蛋白编码基因的测序深度及覆盖度
【直播】我的基因组77：批量计算每个蛋白编码基因的测序深度及覆盖度
生信技能树
发布于 2018-03-08 14:20:56
发布于 2018-03-08 14:20:56
目前我使用的仍然是hg19系统的参考基因组，所以就在gencode数据库里面下载了基于hg19的gtf注释文件，并格式化如下：
head ~/reference/gtf/gencode/protein_coding.hg19.position
我们论坛有专门的教程讲解如何格式化，得到每个基因组的起始终止坐标，就不在此赘述啦（根据gtf格式的基因注释文件得到人所有基因的染色体坐http://www.biotrainee.com/thread-472-1-1.html）
然后我们可以利用大名鼎鼎的bedtools对所有基因批量算出GC含量，如下：
bedtools nuc -fi ~/reference/genome/hg19/hg19.fa -bed ~/reference/gtf/gencode/protein_coding.hg19.position |cut -f 1-4,5,13 >hg19.protein_coding.gc.txt results.txt基因长度只是个附带品咯，因为我本来就有基因的起始终止坐标，所以说长度几乎等于是已知的咯。
bedtools的nuc命令还有给出其它信息，我们并不需要，就取第5,13列即可，基本的shell语法大家需要自己学一点，别看了我的直播这么久，还问那些基础问题。
之前我们讲过samtools的depth用法，很容易就可以根据我们拿到的基因起始终止坐标信息来批量依次提取每个基因的被测序的长度，平均测序深度，还有平均测序深度的方差！
脚本如下：
cat ~/reference/gtf/gencode/protein_coding.hg19.position |while read id
do
arr=($id)
echo ${arr[0]}:${arr[1]}-${arr[2]}
echo -n ${arr[3]} >>results.txt
echo -n -e '\t' >>results.txt
samtools depth -r ${arr[0]}:${arr[1]}-${arr[2]} ../bamFiles/jmzeng.filter.rmdup.bam | \
perl -alne '{$h{$.}=$F[2];$sum+=$F[2]}END{if($.>0){$mean=$sum/$.;$cum+=($_-$mean)**2 foreach values %h;$cum/=$.;print "$.\t$mean\t$cum"}else{print "0\t0\t0"}}' >>results.txt
done
这个脚本很简单，主要是对samtools的depth的输入进行简单的统计而已。
我们可以从统计的结果看到有的基因覆盖度极高，但有的基因覆盖度却很低，这是为什么呢？下一讲我们就简单的解析一下蛋白编码基因的测序深度以及覆盖度吧！
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2017-05-15，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
相关产品与服务
云直播
云直播（Cloud Streaming Services，CSS）为您提供极速、稳定、专业的云端直播处理服务，根据业务的不同直播场景需求，云直播提供了标准直播、快直播、云导播台三种服务，分别针对大规模实时观看、超低延时直播、便捷云端导播的场景，配合腾讯云视立方·直播 SDK，为您提供一站式的音视频直播解决方案。