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

相关文章

来自专栏简书专栏

Scipy入门

标题中的英文首字母大写比较规范,但在python实际使用中均为小写。 建议读者安装anaconda,这个集成开发环境自带了很多包。 作者推荐到2018年8月...

471
来自专栏移动开发面面观

OpenGL ES——打光

935
来自专栏北京马哥教育

Kmeans聚类代码实现及优化

云豆贴心提醒,本文阅读时间6分钟 这篇文章直接给出上次关于Kmeans聚类的篮球远动员数据分析案例,最后介绍Matplotlib包绘图的优化知识。 希望这篇文...

2815
来自专栏大数据挖掘DT机器学习

利用word2vec对关键词进行聚类

按照一般的思路,可以用新闻ID向量来表示某个关键词,这就像广告推荐系统里面用用户访问类别向量来表示用户一样,然后就可以用kmeans的方法进行聚类了。不过对于新...

30810
来自专栏生信技能树

一篇文章学会ChIP-seq分析(下)

写在前面:《一篇文章学会ChIP-seq分析(上)》《一篇文章学会ChIP-seq分析(下)》为生信菜鸟团博客相关文章合集,共九讲内容。带领你从相关文献解读、资...

4617
来自专栏生信宝典

NGS基础 - FASTQ格式解释和质量评估

FASTQ文件格式和命名 高通量测序之后用于下游分析的数据一般存储在FASTQ文件中。为了节省空间,又不影响下游使用,也一般用gzip压缩的格式。 单端测序每个...

2055
来自专栏生信宝典

GSEA富集分析 - 界面操作

GSEA定义 Gene Set Enrichment Analysis (基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋...

2608
来自专栏菩提树下的杨过

“AS3.0高级动画编程”学习:第四章 寻路(AStar/A星/A*)算法 (下)

在前一部分的最后,我们给出了一个寻路的示例,在大多数情况下,运行还算良好,但是有一个小问题,如下图: ? 很明显,障碍物已经把路堵死了,但是小球仍然穿过对角线跑...

1819
来自专栏前端架构

深入理解CSS3 gradient斜向线性渐变

提问,使用CSS3 gradient渐变,在一个400*300的div层上实现一个(100px, 100px)到(200px, 200px)由红到黄的斜向线性渐...

692
来自专栏数据小魔方

R语言学习笔记——柱形图

今天分享R语言中的柱形图,所有图表语法都基于ggplot2包中的ggplot函数完成 。 其实R语言本身就带有各种作图函数,比如plot、bar、pie等,而且...

40613

扫描关注云+社区