专栏首页生信技能树学徒考核-计算wes数据的全部外显子的平均测序深度

学徒考核-计算wes数据的全部外显子的平均测序深度

如果学徒之后跑流程,那其实前途很有限,所以我安排了一个随机任务,考核他们查资料解决问题的能力。我在Published: 04 April 2012 文章, The clonal and mutational evolution spectrum of primary triple-negative breast cancers 看到了一个有趣的图。

首先走wes流程拿到bam文件

这个我们多次讲解了,略,大家自行前往B站看WES视频:

然后根据CCDS数据库拿到人类全部exon的坐标在生信技能树早期教程我也多次讲解过,如何根据CCDS数据库文件,来制作如下BED格式的人类外显子坐标记录文件:

$ head hg38.exon.bed
chr1    69090   70007   OR4F5   0   +
chr1    450739  451677  OR4F29  0   +
chr1    685715  686653  OR4F16  0   +
chr1    801942  802433  LINC00115   0   +
chr1    925941  926012  SAMD11  0   +
chr1    930154  930335  SAMD11  0   +
chr1    931038  931088  SAMD11  0   +
chr1    935771  935895  SAMD11  0   +
chr1    939039  939128  SAMD11  0   +
chr1    939274  939459  SAMD11  0   +

使用samtools工具对exon坐标全部碱基计算覆盖深度

很简单的命令:

~/miniconda2/envs/WES/bin/samtools depth -b hg38.exon.bed a5.sort.bam > /tmp/tmp.depth
$ head tmp.depth
chr1    69091   5
chr1    69092   5
chr1    69093   5
chr1    69094   5
chr1    69095   4
chr1    69096   4
chr1    69097   4
chr1    69098   4
chr1    69099   4
chr1    69100   4

使用bedtools把碱基覆盖深度归属于exon

可以看到每个exon的所以坐标都是有测序深度的,这个文件目前是几千万行!

chr1    69090   70007   OR4F5   0   +   chr1    69091   69091   5
chr1    69090   70007   OR4F5   0   +   chr1    69092   69092   5
chr1    69090   70007   OR4F5   0   +   chr1    69093   69093   5
chr1    69090   70007   OR4F5   0   +   chr1    69094   69094   5
chr1    69090   70007   OR4F5   0   +   chr1    69095   69095   4
chr1    69090   70007   OR4F5   0   +   chr1    69096   69096   4
chr1    69090   70007   OR4F5   0   +   chr1    69097   69097   4
chr1    69090   70007   OR4F5   0   +   chr1    69098   69098   4
chr1    69090   70007   OR4F5   0   +   chr1    69099   69099   4
chr1    69090   70007   OR4F5   0   +   chr1    69100   69100   4

对exon进行汇总

每个坐标的测序深度取平均值即可,可以写一个简短的perl脚本,或者直接读入该文件到R语言,总之对20多万个外显子都计算一个平均测序深度即可。

绘制boxplot

这个是最简单了,参考文献里面的一百多个wes样本合并的boxplot。

课程内容

生信-R语言入门

GEO数据库挖掘

生信-LINUX基础

转录组课题设计和流程分析

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-04

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 对参考基因组按照200k分区间统计测序深度及GC含量

    http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

    生信技能树
  • 你以为的可能不是你以为的

    最近生信技能树管理员小朋友XZG跟我炫耀他植物的简化基因组的gvcf模式,两百个测序数据,我一直没用过这个gvcf功能,因为的确没有需求。癌症研究,关注的主要是...

    生信技能树
  • (14)不同基因坐标转换-生信菜鸟团博客2周年精选文章集

    主流有3个,我只介绍了两个: 用crossmap代替liftover做基因组坐标转换 liftover基因组版本直接的coordinate转换 其实国际三大主流...

    生信技能树
  • 你以为的可能不是你以为的

    最近生信技能树管理员小朋友XZG跟我炫耀他植物的简化基因组的gvcf模式,两百个测序数据,我一直没用过这个gvcf功能,因为的确没有需求。癌症研究,关注的主要是...

    生信技能树
  • (14)不同基因坐标转换-生信菜鸟团博客2周年精选文章集

    主流有3个,我只介绍了两个: 用crossmap代替liftover做基因组坐标转换 liftover基因组版本直接的coordinate转换 其实国际三大主流...

    生信技能树
  • Kubernetes物联网边缘工作组发布白皮书:识别边缘的安全问题

    随着物联网和边缘计算成为构建和部署应用程序的可行选择,越来越多的开发者希望在典型的数据中心部署之外使用Kubernetes和其它云原生技术。这些开发者希望能够使...

    CNCF
  • Salesforce实时聊天工具Live Agent介绍

    Salesforce Live Agent是原生的Salesforce工具,可让企业和网站的用户实时的进行聊天。在这篇文章里我们将讨论为什么需要Live Age...

    臭豆腐
  • html5网页结构布局标签

    对于HTML5来讲,在网页结构上标签定义与使用更加语义化,让搜索引擎以及工程师更加迅速理解当前网页的整个重心所在!

    山河木马
  • 如何扩展一个自定义SOP节点

    平安夜祝大家平平安安,以后的文章关于C++语言方面的内容会多一些,不太理解的话就当一乐子看,了解一下Houdini底层架构知识也是好的。能保证的是文章的内容都是...

    企鹅号小编
  • Lepus搭建企业级数据库全方位监控系统

    Lepus(天兔)数据库企业监控系统是一套由专业DBA针对互联网企业开发的一款专业、强大的企业数据库监控管理系统,企业通过Lepus可以对数据库的实时健康和各种...

    小柒2012

扫码关注云+社区

领取腾讯云代金券