查看原文
其他

【直播】我的基因组78:简单解析一下蛋白编码基因的测序深度及覆盖度

2017-05-18 Jimmy 生信技能树

上一讲中,我们对蛋白的编码基因的测序深度和覆盖度进行了统计,其中有的覆盖度很高,有的覆盖度却又很低,针对这个统计出的测序深度及覆盖度,我们就可以做一些简单的统计及分析。



首先,可以看看覆盖度为10%~100%区间的基因都有多少,并可视化,R代码如下:

  1. hist(dat$coverage,breaks =(0:10)/10)

  2. library(ggplot2)

  3. ggplot(dat,aes(x=coverage))+geom_histogram(binwidth = 0.1 )+

  4. stat_bin(binwidth=0.1, geom="text", aes(label=..count..), vjust=-1.5)+

  5. theme_set(theme_set(theme_bw(base_size=20)))+

  6. theme(text=element_text(face='bold'),

  7. axis.text.x=element_text(angle=30,hjust=1,size =15),

  8. plot.title = element_text(hjust = 0.5) ,

  9. panel.grid = element_blank()

  10. )


很明显,大部分基因(18800/19735=95.26%)都是100%覆盖的,只有少部分基因覆盖度不完整。

值得注意的是,居然有295个基因是完全没有被覆盖到,这个现象值得深究。


我们也顺便看看GC含量跟测序深度的关系。

当然,这里仅仅是选择那些覆盖度大于90%的基因,还有不考虑X,Y,MT等非常染色体基因咯。

  1. dat_new=subset(dat,coverage>0.9)

  2. dat_new=subset(dat_new,chr %in% paste0('chr',1:22))

  3. plot(dat_new$depth,dat_new$gc)


这个趋势实在是太明显了,GC含量越高,测序深度越低,说明二代测序的硬伤是存在的。因为GC含量高的区域很难PCR扩增。

上图我截掉了测序深度超过100的那些基因,单独显示如下:

这些基因为什么测序深度这么高呢?我的WGS数据全基因组平均测序深度只有45X。

接着回过头看看那230个完全没有被覆盖到的基因吧~

  1. dat_new=subset(dat,coverage ==0)

  2. sort(table(dat_new$chr))

  3. barplot(sort(table(dat_new$chr)))


我看了一下,6号染色体就占了一多半了,很有可能是6号染色体的注释不够完全,而不是我的基因组的问题。因为性染色体就排在后面,它们上面的基因没办法覆盖到这很正常了。

我仔细检查了6号染色体的这些基因,发现很多是orf系列基因,我在我们生信技能树论坛里面曾经发帖提到过这件事情。然后就是一堆主要组织相容性的复合物,这个没什么好说的了,免疫的相关基因本来就乱乱的。还有几个锌指蛋白,几个嗅觉受体蛋白编码基因,还有一些多肽,反正我也不认识。也说不出来个所以然来。


至于性染色体中,X主要是几个cancer/testis antigen family基因,还有cancer/testis antigen family chromosome X open reading frame基因,SPANX家族,X 抗体家族,G抗体家族,还有热激蛋白。Y染色体上面没有被覆盖到的基因,我貌似都不认识呀。

而1号染色体上面覆盖度为0的都是histone cluster基因,为什么它们无法被测序呢?GC含量比较低38%,可是GC含量低,应该是测序深度高呀!基因长度比较短,貌似也不是理由。

4号染色体我检查了一下,都是ubiquitin specific peptidase 17-like family member系列基因,这个很容易理解了,本身这个家族基因注释就很烂,它们内部相似性太高了,比对的时候压根就没办法把reads唯一定位到家族的某个具体基因。所以导致家族某些基因超高深度测序结果,而有些基因根本就没有测序结果。

值得关注的现象应该还有很多,我就不一一解读啦,希望各位读者朋友们可以跟我来信讨论交流感兴趣的地方。

 推荐阅读:

直播我的基因组分析-目录-1-73集

【直播】我的基因组74:快速给测序reads比对到物种

【直播】我的基因组75:为什么我的血液里面有这些菌?

【直播】我的基因组76:用krona对血液全基因组的菌比例可视化

【直播】我的基因组77:批量计算每个蛋白编码基因的测序深度及覆盖度


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

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