查看原文
其他

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

2017-09-11 曾健明 生信技能树

一般情况下,我们得到了测序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. 1    226 exon:43

  2. 2    235 exon:43

  3. 3    246 exon:43

  4. 4    254 exon:43

  5. 5    258 exon:43

  6. 6    262 exon:43

  7. 7    277 exon:43

  8. 8    286 exon:43

  9. 9    298 exon:43

  10. 10    319 exon:43

  • 第一列是该外显子的起始终止坐标,从1到该外显子的长度。每切换一个外显子,坐标从1开始记录。一般的外显子长度在200bp左右,也有一些超长的外显子5kb及以上长度。

  • 第二列是该外显子在该坐标的测序深度,通过samtools的depth命令得到

  • 最后一列是该外显子的标记,从exon:79一直倒推到exon:1,因为该基因在染色体的负链,所以外显子顺序是反着的!这一列信息用来把外显子在ggplot上面分面!

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

  1. if("ggplot2" %in% rownames(installed.packages()) == FALSE) {install.packages("ggplot2")}

  2. library(ggplot2)

  3. args <- commandArgs(trailingOnly = TRUE)

  4. file=args[1]

  5. outpng=sub(".txt",".png",file)

  6. dat=read.table(file)

  7. names(dat)=c('pos','depth','exon')

  8. new_order=paste0('exon:',1:length(unique(dat[,3])))

  9. dat[,3] <- factor(dat[,3] ,new_order )

  10. png(outpng,res=120,width = 1080, height = 1080)

  11. p=ggplot(data=dat,aes(x=pos,y=depth,color=exon))+geom_line()

  12. p=p+facet_wrap(~exon,scales="free_x")

  13. p=p+theme(legend.position='none')

  14. print(p)

  15. dev.off()

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

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

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

  1. bam=$1

  2. file=$(basename $bam )

  3. sample=${file%%.*}

  4. echo $sample

  5. mkdir -p  exon_png_$sample

  6. genes=(DMD TP53 BRCA1 BRCA2 ALK ROS1 EGFR MET BRAF FGFR1 RET KRAS)

  7. for gene in ${genes[@]}

  8. do

  9. perl ~/biosoft/exoncoverage/bin/cal.pl  -r ~/biosoft/exoncoverage/db/hg19_refGene.txt \

  10. -b $bam  -S $gene

  11. done

  12. mv *.png  exon_png_$sample

  13. rm *.txt

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

可以看到我的WGS测序深度在50~60之间,这很正常!因为我就花了这些钱而已。

不过BRCA1的第2,8外显子几乎没有reads覆盖,这个很诡异,但是因为这两个外显子长度小于100bp,很有可能是NGS测序技术的限制,所以我倒是没那么担心。

如果大家有人类的WGS测序数据的,可以用我的方法做一个类似的图,跟我交流一下,我的邮箱是 jmzeng1314@163.com  希望你可以给我一点信心,我很正常。

猜你喜欢

基因组 游记 | 工作资讯 

学习课程 | 好书分享 


菜鸟入门

Linux | Perl | R语言 | 可视化 

R包 | perl模块 | python模块


数据分析

ChIP-seq(上)ChIP-seq(下)RNA-seq | miRNA

WGS,WES,RNA-seq组与ChIP-seq之间的异同


编程实践

第0题 | 探索人类基因组序列 | 最后一题


直播基因组分析

我的基因组 | 解惑帖

一个标准的基因检测报告目录

生信技能树

直播我的基因组分析-目录-持续更新

【直播】我的基因组81:看看我的vcf文件的vaf分布情况

【直播】我的基因组82:如何对maf格式的突变文件统计vaf

批量IGV截图【直播】我的基因组83

从WGS测序得到的VCF文件里面提取位于外显子区域的【直播】我的基因组84

安装snpEFF工具并对VCF文件进行注释【直播】我的基因组85

基因组重测序的unmapped reads assembly探究 【直播】我的基因组86


点击下面的阅读原文可以拿到本次分析的所有代码!

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存