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

相关文章

来自专栏SeanCheney的专栏

《Python数据分析》2nd

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

3768
来自专栏Y大宽

Enrichment Map User guide用户指南

http://www.baderlab.org/Software/EnrichmentMap/UserManual#rnk

2383
来自专栏点滴积累

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

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

3475
来自专栏数据和云

拨开迷雾:关于resize datafile理解的错误

黄廷忠(网名:认真就输) 云和恩墨技术专家 个人博客:http://www.htz.pw/ resize数据文件的时候,常常会报ORA-03297错误,也许很多...

3304
来自专栏生信宝典

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

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

5653
来自专栏高爽的专栏

登录之验证码

       产生验证码,MakeCertPic.java: public class MakeCertPic { // 验证码图片中可以出现的字符集,可...

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

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

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

3546
来自专栏LET

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

3034
来自专栏生信技能树

bedtools 用法大全(一文就够吧)

bedtools等工具号称是可以代替普通的生物信息学数据处理工程师的!我这里用一个专题来讲解它的用法,其实它能实现的需求,我们写脚本都是可以做的,而且我强烈建议...

1.1K6
来自专栏HansBug's Lab

算法模板——Dinic最小费用最大流

实现功能:输入M,N,S,T;接下来M行输入M条弧的信息(包括起点,终点,流量,单位费用);实现功能是求出以S为源点,T为汇点的网络最大流的最小费用 其实相当的...

3556

扫码关注云+社区