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

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

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

hist(dat$coverage,breaks =(0:10)/10)
library(ggplot2)
ggplot(dat,aes(x=coverage))+geom_histogram(binwidth = 0.1 )+
stat_bin(binwidth=0.1, geom="text", aes(label=..count..), vjust=-1.5)+
theme_set(theme_set(theme_bw(base_size=20)))+
theme(text=element_text(face='bold'),
axis.text.x=element_text(angle=30,hjust=1,size =15),
plot.title = element_text(hjust = 0.5) ,
panel.grid = element_blank()
)

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

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

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

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

dat_new=subset(dat,coverage>0.9)
dat_new=subset(dat_new,chr %in% paste0('chr',1:22))
plot(dat_new$depth,dat_new$gc)

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

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

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

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

dat_new=subset(dat,coverage ==0)
sort(table(dat_new$chr))
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唯一定位到家族的某个具体基因。所以导致家族某些基因超高深度测序结果,而有些基因根本就没有测序结果。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信宝典

癌症组织特异性基因怎么找?这是个不错的开始

组织特异性基因(Tissue-specific Genes)是指在不同类型的细胞中特异性表达的基因,其调节细胞特异的形态结构或生理功能。组织特异性基因的表达是理...

142
来自专栏FreeBuf

浅析FPGA现场可编程门阵列

FPGA(Field-Programmable Gate Array),即现场可编程门阵列,它是在PAL、GAL、CPLD等可编程器件的基础上进一步发展的产物。...

1958
来自专栏绿巨人专栏

读书笔记: 博弈论导论 - 16 - 不完整信息的动态博弈 信号传递博弈

2587
来自专栏程序员互动联盟

【编程心里】编程大牛教你正确的学习心态

小明问大师,大师我已经开始学习c语言编程了为什么感觉我只会用他做数学题,而不能写自己想写的游戏呢? 大师看着地上的教学文章不说话; 小明说大师你是让我静心学习之...

3535
来自专栏大数据钻研

2017年前端工程师应该学习什么

在我们所生活的这个快节奏的世界里,人们都倾向于把自己的时间用在进行一些新的创造上,然后再互联网上讨论它们。 我并不是说不该这样做,而是我认为我们应该适当的慢一...

3296
来自专栏ATYUN订阅号

IBM团队开发新的AI算法,可以过滤侮辱性语言并以礼貌用语来代替

这是记者们经常重复的一句话:从不读评论。评论部分的内容,可能会是互联网上最黑暗的地方之一,那些毫无根据的侮辱和尖锐的批评像混战中的子弹一样。

1014
来自专栏PPV课数据科学社区

当禅师遇到一位理科生,后来禅师疯了!!知识无极限!!

老禅师忙着给各种年轻人指导人生,总是用一些模凌两可的语句,想着想着你自己似乎就想透了。但当满口心灵鸡汤的老禅师遇上理科生……于是有了下面的对话↓↓↓ 1、青...

3333
来自专栏生信技能树

【直播】我的基因组80:为什么有些基因的内部测序深度差异如此大

这一讲里,我们依旧根据统计的基因测序的深度进行一下讨论,来看看为什么有些基因的内部测序深度差异如此大? 在前面我们的计算中,s列表示的是基因的每一个坐标的...

3277
来自专栏PaddlePaddle

【AI核心技术】课程二十一:神经图灵机—控制器

UAI与PaddlePaddle联合推出的【AI核心技术掌握】系列课程持续更新中!

563
来自专栏企鹅号快讯

如何在短时间迅速提升Python功力?Python应该怎样去学习呢?

最近有一段时间没有看后台的留言,。昨天我打开后台一看有很多同学给我留言,其中有5条是问我关于如何快一点提高Python功力的相关问题~~ 确实当你学了Pytho...

1959

扫描关注云+社区