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

相关文章

来自专栏Python小屋

Python+tkinter实现任意多层级关系的组合框

正好自己要用,就想了个简单思路实现了一下,假设现在需要一个组合框来显示多个层级关系的单位名称供用户选择,如果只有固定的两级关系当然可以使用两个组合框联动来实现,...

3175
来自专栏生信宝典

GO、GSEA富集分析一网打进

富集分析是生物信息分析中快速了解目标基因或目标区域功能倾向性的最重要方法之一。其中代表性的计算方式有两种: 一是基于筛选的差异基因,采用超几何检验判断上调或下调...

28610
来自专栏转载gongluck的CSDN博客

FFmpeg菜鸡互啄#第4篇#音频解码

解码过程 音频解码跟上一篇的视频解码过程是一样的:打开输入文件,查找音频流,打开解码器,循环读帧解码帧,关闭解码器,关闭输入文件。 Code #define _...

2946
来自专栏何俊林

Android Multimedia框架总结(二十)MediaCodec状态图及Codec与输入/输出Buffer过程(附实例)

前言:前面几节都是介绍Camera2相关,对于Camera2预览把图像显示在SurfaceView上,还有录像时,时时刷新当前图像区域。追溯到最早介绍的Medi...

2168
来自专栏偏前端工程师的驿站

动手写个数字输入框1:input[type=number]的遗憾

前言  最近在用Polymer封装纯数字的输入框,开发过程中发现不少坑,也有很多值得研究的地方。本系列打算分4篇来叙述这段可歌可泣的踩坑经历: 《动手写个数字输...

2115
来自专栏大闲人柴毛毛

高质量编程的金玉良言——依赖倒转原则

生活中的例子: 电脑的品牌有很多,但电脑中的所有部件都有标准的接口,不同的厂家都是按照标准去生产各个部件,这些部件的内部实现不同,但接口都是一样的,这样的话,如...

2667
来自专栏何俊林

Android Multimedia框架总结(九)Stagefright框架之数据处理及到OMXCodec过程

不知不觉到第九篇了,感觉还有好多好多没有写,路漫漫其修远兮 ,吾将上下而求索。先说福利吧,此前在关于我, ? 曾说过,不定期搞活动,vip,书啥的,都可以有,...

1896
来自专栏鸿的学习笔记

如何设计一个良好的流系统?(下)

在Streaming 101中,作者引入了窗口和时间的概念,在本文中,作者为了解决流处理系统无法精确的处理结果的问题,提出了下面三个概念:

231
来自专栏落影的专栏

使用VideoToolbox硬解码H.264

前言 使用VideoToolbox硬编码H.264 在上一篇的硬编码简单介绍了H.264和VideoToolbox以及如何使用VideoToolbox硬编码从...

4225
来自专栏小樱的经验随笔

【批处理学习笔记】第二十四课:直接传递

    直接传递参数,即在使用call命令时,不使用任何参数,在子函数或子批处理里面直接对主函数(也称父批处理)里面的变量进行修改。这跟汇编语言里面的参数传递方...

2383

扫描关注云+社区