探究某个基因的外显子覆盖度情况【直播】我的基因组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的转录本的其中一个外显子看起来是缺失了。如下：