探究某个基因的外显子覆盖度情况【直播】我的基因组87

一般情况下,我们得到了测序reads在基因组的比对情况文件bam格式的,里面的信息非常多,如果我想特定的查看某个基因的情况,那么我们可以选择IGV等可视化工具,但它并不是万能的,因为即使是一个基因,它也会有多个转录本,多个外显子。以前我写过批量IGV截图(点击直达),但是大部分基因的长度都超过了37Kb,超过了IGV的窗口浏览限制。而且我们也不需要知道该基因上面比对成功的所有reads信息,太复杂了,我们只需要知道基因上面各个部位的测序深度即可,而且基因上面比较重要的就是外显子了,被一个个内含子隔离开来。

所以,我们可以画它的外显子覆盖图,下面是一个例子:横坐标是外显子的长度,纵坐标是测序深度,每一个小图都是一个外显子

DMD外显子覆盖度情况

根据这个图,我们就可以很明显的看出,DMD基因NM_000109转录本的1,10-17号外显子缺失,用IGV一个个的看这些外显子区域,是同样的结果!可能是芯片捕获不到,也可能是样本本身变异,造成的大片段缺失。但是这个图的信息就非常有用!(当然,这个肯定不是我的啦,我很正常的哦)

PS:请忽略上图的外显子不是按照数值的大小排序,只是因为当初我对ggplot还不是很熟悉,不知道调整factor就可以改变出图的顺序。

那么,我们该如何画这样的图呢?

首先,我们需要找到需要探究的基因的全部转录信息,及外显子信息!

这里的hg19_refGene.txt我直接从annovar的数据目录拿到的。可以看到一个基因有多个转录本,每个转录本的外显子个数不一样。如下:

那么,我们根据这个信息,就可以判断该基因的每个外显子的起始终止位点啦!

然后用samtools的depth命令去找这个基因的全部片段的测序深度信息

最后再格式化成下面的三列数据

1    226 exon:432    235 exon:433    246 exon:434    254 exon:435    258 exon:436    262 exon:437    277 exon:438    286 exon:439    298 exon:4310    319 exon:43
  • 第一列是该外显子的起始终止坐标,从1到该外显子的长度。每切换一个外显子,坐标从1开始记录。一般的外显子长度在200bp左右,也有一些超长的外显子5kb及以上长度。
  • 第二列是该外显子在该坐标的测序深度,通过samtools的depth命令得到
  • 最后一列是该外显子的标记,从exon:79一直倒推到exon:1,因为该基因在染色体的负链,所以外显子顺序是反着的!这一列信息用来把外显子在ggplot上面分面!

最后根据这个txt文档,用R语言,很容易就画出上面那样的图片了!R语言程序是:

if("ggplot2" %in% rownames(installed.packages()) == FALSE) {install.packages("ggplot2")}library(ggplot2)args <- commandArgs(trailingOnly = TRUE)file=args[1]outpng=sub(".txt",".png",file)dat=read.table(file)names(dat)=c('pos','depth','exon')new_order=paste0('exon:',1:length(unique(dat[,3])))dat[,3] <- factor(dat[,3] ,new_order )png(outpng,res=120,width = 1080, height = 1080)p=ggplot(data=dat,aes(x=pos,y=depth,color=exon))+geom_line()p=p+facet_wrap(~exon,scales="free_x")p=p+theme(legend.position='none')print(p)dev.off()

这个是修正之后的R代码,外显子的顺序ok了。

比如下面这个捕获测序的,可以很明显的看到外显子区域测序深度超级高!但是第一个和最后一个外显子很明显测序深度不足,有可能是设计这个panel的人认为第一个和最后一个外显子是UTR区域,不予捕获。

用这个代码画一下我的WGS数据

bam=$1file=$(basename $bam )sample=${file%%.*}echo $samplemkdir -p  exon_png_$samplegenes=(DMD TP53 BRCA1 BRCA2 ALK ROS1 EGFR MET BRAF FGFR1 RET KRAS)for gene in ${genes[@]}doperl ~/biosoft/exoncoverage/bin/cal.pl  -r ~/biosoft/exoncoverage/db/hg19_refGene.txt \-b $bam  -S $genedonemv *.png  exon_png_$samplerm *.txt

可以看到,我只是挑选了部分基因来画图,大部分都挺好的, 全部覆盖到了,但是有一个BRCA1的转录本的其中一个外显子看起来是缺失了。如下:

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

原文发表时间:2017-09-11

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏懒人开发

(10.6)James Stewart Calculus 5th Edition:Conic Sections in Polar Coordinates

952
来自专栏深度学习入门与实践

【原】Spark之机器学习(Python版)(一)——聚类

kmeans聚类相信大家都已经很熟悉了。在Python里我们用kmeans通常调用Sklearn包(当然自己写也很简单)。那么在Spark里能不能也直接使用s...

88010
来自专栏简书专栏

基于tensorflow、CNN、清华数据集THUCNews的新浪新闻文本分类

tensorflow是谷歌google的深度学习框架,tensor中文叫做张量,flow叫做流。 CNN是convolutional neural netwo...

4731
来自专栏C#

开源免费的.NET图像即时处理的组件ImageProcessor

   承接以前的组件系列,这个组件系列旨在介绍.NET相关的组件,让大家可以在项目中有一个更好的选择组件的介绍绝对不是一篇文章可以叙述完的,因为一个组件是经过开...

2548
来自专栏量化投资与机器学习

【深度学习量化投资】RNNs在股票价格预测的应用基于Keras

前言 RNN和LSTMs在时态数据上表现特别好,这就是为什么他们在语音识别上是有效的。我们通过前25天的开高收低价格,去预测下一时刻的前收盘价。每个时间序列通过...

1.3K6
来自专栏机器学习从入门到成神

Python3读取深度学习CIFAR-10数据集出现的若干问题解决

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/sinat_35512245/articl...

862
来自专栏腾讯云TStack专栏

《 大话 Ceph 》 之 CRUSH 那点事儿

《大话 Ceph 》系列文章通过通俗易懂的语言并结合基础实验,用最简单的描述来讲解 Ceph 中的重要概念。让读者对分布式存储系统有一个清晰的理解。

1.9K5
来自专栏贾志刚-OpenCV学堂

手撕OpenCV源码之filter2D(一)

在上篇的GaussianBlur中提到,gaussianBlur使用的是filter2D的实现,因此上篇仅仅描述了高斯滤波器的生成细节,并没有针对滤波的计算细节...

1601
来自专栏生信宝典

Bedtools使用简介

3614
来自专栏生信技能树

【直播】我的基因组73:在IGV看看indel是啥样子

前面我们特意用scalpel软件来找indel,期待它会有一些出彩的表现,当然我还没来得及比较它找到的INDEL跟GATK等工具区别在哪里,不过我们先在IGV里...

4139

扫码关注云+社区