【直播】我的基因组77:批量计算每个蛋白编码基因的测序深度及覆盖度

目前我使用的仍然是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的输入进行简单的统计而已

我们可以从统计的结果看到有的基因覆盖度极高,但有的基因覆盖度却很低,这是为什么呢?下一讲我们就简单的解析一下蛋白编码基因的测序深度以及覆盖度吧!

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-15

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信技能树

【直播】我的基因组65:看看哪些基因的突变较多,哪些较少

全基因组分析后的vcf突变文件记录了四百多万个位点,前面我们讲到了如何把它们注释到dbSNP数据库ID,一般来说有注释的位点也就顺便注释到了基因,所以可以简单写...

3119
来自专栏前端大白专栏

基于mpvue开发微信小程序(项目已开源)

1816
来自专栏星汉技术

Scala简介和安装

2676
来自专栏Python小屋

Python使用matplotlib进行可视化时精确控制图例位置

在进行数据可视化或者科学计算可视化时,显示图例会显得很高大上,但是如果能够精确控制图例的显示位置,无疑会显得档次更高。 本文以matplotlib.pyplot...

2616
来自专栏生信技能树

OMIM使用简要说明【论坛精选优秀帖】

.OMIM 为“0nline MendelianInheritance in Man”的简称,它通过对新的病症分类并命名、收录表型和相关病因基因的关系来收录人类...

32910
来自专栏美团技术团队

Android热更新方案Robust

美团•大众点评是中国最大的O2O交易平台,目前已拥有近6亿用户,合作各类商户达432万,订单峰值突破1150万单。美团App是平台主要的入口之一,O2O交易场景...

3719
来自专栏前端大白专栏

基于mpvue开发微信小程序(项目已开源)

花了两周时间,我的微信小程序终于开发完了(平时上班,基本上都是业余时间开发的). 下面来介绍一下项目的功能以及结构. 用到的技术栈 vue2+weui+es6;...

4769
来自专栏窗户

shell编程/字库裁剪(1)——想法

  我写这个帖子的意图,在于三个:   1.用代码生成代码的思维。   2.shell编程的思路。   3.裁剪字库的具体程序。   我打算分为三节来说:   ...

2019
来自专栏程序员互动联盟

程序语言变形记

随着科技的发展我们生活中接触到的应用程序越来越多,它给我们的生活带来了很大的便利。移动端安卓,苹果大肆横行;pc上QQ,浏览器大行天下。我们在享受这些软件给我们...

3095
来自专栏刘望舒

热修复原理之热修复框架对比和代码修复

题外话 今天听到了著名物理学家史蒂夫霍金去世消息,潸然泪下。作为一个天文迷,感谢霍金带给我的那些天文知识。十分敬佩霍金的身残志坚,他在全身瘫痪无法言语情况下仍旧...

3084

扫码关注云+社区