【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏利炳根的专栏

学习笔记TF063:TensorFlow Debugger

TensorFlow Debugger(tfdbg),TensorFlow专用调试器。用断点、计算机图形化展现实时数据流,可视化运行TensorFlow图形内部...

65500
来自专栏LET

谈谈3D Tiles(2):数据结构

44240
来自专栏数据结构与算法

51NOD 1185 威佐夫游戏 V2(威佐夫博弈)

1185 威佐夫游戏 V2 基准时间限制:1 秒 空间限制:131072 KB 分值: 0 难度:基础题 有2堆石子。A B两个人轮流拿,A先拿。每次可以从一...

37060
来自专栏和蔼的张星的图像处理专栏

DSP图像处理

最近着手把CSK移植到DSP中,先看一些DSP中图像处理的一些例子,第一件事当然就是怎么把图像数据倒入CCS工程中了,去年倒是用过一点CCS,再拿起来已经忘得差...

28520
来自专栏生信宝典

分子对接简明教程 (4)

文件格式解释 PDB文件 (详细格式描述) 基本信息部分 HEADER记录: 包括分子的分类、提交日期、PDB ID TITLE记录: 为该结构的描述,如果有多...

45270
来自专栏吉浦迅科技

DAY46:阅读Surface Reference API

reads the CUDA array bound to the one-dimensional surface reference surfRef usin...

11850
来自专栏SeanCheney的专栏

《Python数据分析》2nd

《Python数据分析》(Python for Data Analysis, 2nd Edition)第二版出了,目前还没有中文版,这版的代码适用于Python...

41580
来自专栏点滴积累

geotrellis使用(二十八)栅格数据色彩渲染(多波段真彩色)

目录 前言 实现过程 总结 一、前言        上一篇文章介绍了如何使用Geotrellis渲染单波段的栅格数据,已然很是头疼,这几天不懈努力之后工作又进了...

37550
来自专栏Java架构沉思录

聊聊设计模式之策略模式

前言 这几天大部分同学应该都过完年陆陆续续回到工作岗位了,说到过年,最开心的莫过于与家人团聚了,当然除了与家人团聚,最令人振奋的事情就是发年终奖了。说到发年终...

33270
来自专栏生信宝典

如何获取目标基因的转录因子(上)——Biomart下载基因和motif位置信息

科研过程中我们经常会使用Ensembl(http://asia.ensembl.org/index.html) 网站来获取物种的参考基因组,其中BioMart工...

1.1K30

扫码关注云+社区

领取腾讯云代金券