【直播】我的基因组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唯一定位到家族的某个具体基因。所以导致家族某些基因超高深度测序结果，而有些基因根本就没有测序结果。