探究某个基因的外显子覆盖度情况

2017-09-11 08:00:36生信技能树

一般情况下，我们得到了测序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:43
2    235 exon:43
3    246 exon:43
4    254 exon:43
5    258 exon:43
6    262 exon:43
7    277 exon:43
8    286 exon:43
9    298 exon:43
10    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=$1
file=$(basename $bam )
sample=${file%%.*}
echo $sample
mkdir -p  exon_png_$sample
genes=(DMD TP53 BRCA1 BRCA2 ALK ROS1 EGFR MET BRAF FGFR1 RET KRAS)
for gene in ${genes[@]}
do
perl ~/biosoft/exoncoverage/bin/cal.pl  -r ~/biosoft/exoncoverage/db/hg19_refGene.txt \
-b $bam  -S $gene
done
mv *.png  exon_png_$sample
rm *.txt