探究某个基因的外显子覆盖度情况【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏hightopo

根据矩阵变化实现基于 HTML5 的 WebGL 3D 自动布局

873
来自专栏AI科技评论

开发 | 如何利用微信监管你的TF训练

AI科技评论按:本文作者Coldwings,AI科技评论获其授权发布。 之前回答问题【在机器学习模型的训练期间,大概几十分钟到几小时不等,大家都会在等实验的时候...

2978
来自专栏数据小魔方

时间管理的工具——甘特图(Gantt chart)

今天跟大家分享一种用作时间管理的工具——甘特图(Gantt Chart)。 ▽▼▽ 这种图表的制作理念非常简单,就是通过设定项目开始时间和持续时间,利用堆积条形...

3767
来自专栏破晓之歌

通过 Python 接口使用 OpenCV 原

1.通过 Python 接口使用 OpenCV :https://blog.csdn.net/mzpmzk/article/details/68952160

592
来自专栏AI研习社

如何利用微信监管你的TF训练?

之前回答问题【在机器学习模型的训练期间,大概几十分钟到几小时不等,大家都会在等实验的时候做什么?(http://t.cn/Rl8119m)】的时候,说到可以用微...

3424
来自专栏埋名

CSS样式clip-path绘制图形裁剪图片

发现clip-path绘制图形的网站 http://bennettfeely.com/clippy/ 专门介绍clip-path的使用,几乎任何形状都可以做出来...

422
来自专栏HT

根据矩阵变化实现基于 HTML5 的 WebGL 3D 自动布局

在数学中,矩阵是以行和列排列的数字,符号或表达式的矩形阵列,任何矩阵都可以通过相关字段的标量乘以元素。矩阵的主要应用是表示线性变换,即f(x)= 4 x等线性函...

1985
来自专栏進无尽的文章

UI篇-Layer几个关键点补充

强大的UIView是基于 CALayer实现的,它的重要性不言而喻,相信大家也都有自己的研究和理解,今天这片文章里的内容是几个关键点的补充。

551
来自专栏杨建荣的学习笔记

在oracle实践学习位运算 第一篇(r4笔记第55天)

今天无意中看到了谭浩强先生的<<c程序设计>>这本书,虽然c语言都是很多年前学过的东西了,但是看起来亲切,实际用起来陌生,很多的概念都已经很模糊了,记得上大学时...

2624
来自专栏生信技能树

第3篇:用MACS2软件call peaks

Peak calling即利用计算的方法找出ChIP-seq或ATAC-seq中reads富集的基因组区域。

1634

扫码关注云+社区