【直播】我的基因组50:从测序深度和位点间距来看SNV分布情况

今天的我们,还是继续探究那一个困扰我这么久的问题。为什么我作为堂堂正正的男性,明明X,Y染色体都只有一条,可是却测到了那么多的杂合突变的问题。

在之前,我们在QC阶段详细的探究了X,Y染色体的覆盖度和测序深度,其中X的平均测序深度才16x,而Y却高达60x,我们完全有理由怀疑测序深度对SNV的准确性影响甚大!而且Y染色体总共长度才60M,就有一半是N碱基,有效长度就30M不到,却找到了近3万个SNV,这有着很明显的问题,太密集了~

所以从测序深度和位点间距来看SNV分布情况是非常有必要的!

对于我的X染色体来说,纯合SNV(91588个)占绝大多数,所以问题不大,因为我只有一条X染色体嘛!我重点就分析一下那18424个杂合的SNV看看是什么情况。

我把VCF文件根据纯合和杂合的SNV分开统计了测序深度,并可视化如下:

ChrX杂合SNV深度统计(横坐标是测序深度,纵坐标是SNV个数)

ChrX纯合SNV深度统计(横坐标是测序深度,纵坐标是SNV个数)

很明显,X染色体的纯和SNV是正常的,所以测序深度也非常正常,集中在10~30,超过40的基本就没有多少了。但是对于不正常的那些杂合SNV来说,就很可怕了,测序深度很明显平均高于纯合SNV。而且也高于X染色体的平均测序深度(16X)

那么我们再看看Y染色体,统计并且可视化如下:

ChrY杂合SNV深度统计(横坐标是测序深度,纵坐标是SNV个数)

ChrY纯合SNV深度统计(横坐标是测序深度,纵坐标是SNV个数

同样,纯和SNV是正常的,所以测序深度也非常正常,集中在10~30x,但是超过100的居然还有一大堆!(这些位点太可疑了)而且对于那些不正常的杂合位点来说,很明显测序深度都远高于40x。

接下来我还探究了X,Y染色体的SNV的间距!

这个图对很多人来说比较难以理解,如果需要自己动手实现,要仔细研究我的R代码,其实就是把SNV的坐标提取出来,按照大小排序,然后相邻的坐标之间是有距离的,那么把这些间距拿出来就可以画一个箱线图了,如果箱线图都压缩在一起,就说明大部分SNV的间距实在是太小了,比如X染色体的杂合SNV,简直是一个挨着一个的,完全不符合常理。

也很明显可以看到,我的性染色体的杂合snp位点距离太近了(尤其是X染色体的杂合位点,简直全部凑在了一起,看不出箱线图了),相比纯合来说!

从IGV也可以看到这样的现象大量存在,我随意展示一个一个基因的一个片段reads覆盖截图:

这个基因就这么一个小片段,上面全部是杂合的SNV!

很明显,我可以确定这些杂合SNV是因为人类的X,Y染色体片段大量同源,这样测序得到的reads,必然可以在两个染色体都找到定位,我们的测序策略是PE150,这个长度不足以把片段定位到具体某一条性染色体上面。对于那些multiple mapping的情况,BWA的选择策略是随机给予一个定位即可,那样就导致X,Y染色体片段大量同源变成了2倍体,而不只是理论上的单倍体啦!

我也找到了理论支持,这样的区域叫做PAR(pseudoautosomal region ),包括 PAR1, PAR2,PAR3( 维基百科有详细介绍)

虽然我们对比对的bam文件进行了基本的过滤,比如PCR的duplication,multiple mapping情况,低质量的比对。不过我不太相信,宾夕法尼亚州立大学教授Istvan Albert,他在https://www.biostars.org/p/76893/ 说到按照quality来filter,因为BWA会给多个定位的比对情况质量值为0。虽然软件说明书也是这样写的,但事实上,这样过滤之后,还是有很多reads含有XA:Z,而这个tag意味着有着其它比对选择。

所以multiple mapping的quality必然为0这个预设是值得怀疑的。

Y染色体上面的unique mapping reads共有2.55m, 而multiple mapping reads有1.197m,X染色体上面的unique mapping reads共有16.6m, 而multiple mapping reads有0.971m(m是百万条的意思).

也就是说X,Y染色体上面仍然各有百万条reads是可以多比对情况的,这些reads所定位的区域找到的SNV,都是不可靠的!正是因为这些多比对情况的reads瞎定位,导致了X,Y染色体的测序深度差异如此之大~

上面的统计shell代码是:

bcftools  view -r chrX jmzeng.bcftools.vcf.gzbcftools  view -r chrY jmzeng.bcftools.vcf.gzcat chrX.bcftools.vcf |grep -v "^#" |grep -v INDEL |perl -alne '{/DP=(\d+);/;$depth=$1;@tmp=split(":",$F[9]);print "$F[1]\t$depth\t$tmp[0]"}' >chrx.txtcat chrY.bcftools.vcf |grep -v "^#" |grep -v INDEL |perl -alne '{/DP=(\d+);/;$depth=$1;@tmp=split(":",$F[9]);print "$F[1]\t$depth\t$tmp[0]"}' >chry.txt

对输出的chrx.txt和chry.txt文件进行可视化的R代码是:

bar_plot <- function(a,prefix){  png(paste0(prefix,'.png'),width = 800)  library(calibrate)  tmp = hist(a,breaks = c(1, 10*(1:10),max(a)))  plot(tmp$counts, xaxt='n',type = 'h',lwd=15,xlab = 'depth-region',       ylab='counts',ylim = c(0,1.1*max(tmp$counts)),xlim=c(1,11.5),main = prefix       );  axis(side=1,        labels=paste(c(1, 10*(1:10)),c(10*(1:10),max(a)),sep = '~'),       at=1:11                         )   textxy(1:11,tmp$counts+0.1*max(tmp$counts),tmp$counts,cex = 2,offset=0.3)  dev.off()}chrX=read.table('chrx.txt')colnames(chrX)=c('pos','depth','type')table(chrX$type)# chrX <- chrX[chrX$depth>10,]# table(chrX$type)bar_plot(chrX[chrX$type=='1/1',2],"chrX_hom")bar_plot(chrX[chrX$type=='0/1',2],"chrX_het")chrY=read.table('chrY.txt')colnames(chrY)=c('pos','depth','type')table(chrY$type)# chrY <- chrY[chrY$depth>10,]# table(chrY$type)bar_plot(chrY[chrY$type=='1/1',2],"chrY_hom")bar_plot(chrY[chrY$type=='0/1',2],"chrY_het")par(mfrow=c(2,2))a = chrX[chrX$type=='1/1','pos'] a_diff = diff(a)a_diff = a_diff[ a_diff <10000]boxplot(a_diff,main='chrX-hom')a = chrX[chrX$type=='0/1','pos'] a_diff = diff(a)a_diff = a_diff[ a_diff <10000]boxplot(a_diff,main='chrX-het')a = chrX[chrY$type=='1/1','pos'] a_diff = diff(a)a_diff = a_diff[ a_diff <10000]boxplot(a_diff,main='chrY-hom')a = chrX[chrY$type=='0/1','pos'] a_diff = diff(a)a_diff = a_diff[ a_diff <10000]boxplot(a_diff,main='chrY-het')

文:Jimmy

图文编辑:吃瓜群众

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

原文发表时间:2017-01-24

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏量化投资与机器学习

【顶级资源】掌握线性代数为机器学习打下坚实基础!

线性代数是数学领域,也是机器学习领域重要的支柱。对于初学者来说,要想学好机器学习,线性代数的掌握是必不可少的,也可以说是十分十分重要的。 春节后的第一天,公众号...

4387
来自专栏生信小驿站

文献翻译:Statistical Approaches for Gene Selection, Hub Gene Identification and Module Interaction in...

信息基因的选择是基因表达研究中的重要问题。基因表达数据的小样本量和大量基因特性使选择过程复杂化。此外,所选择的信息基因可以作为基因共表达网络分析的重要输入。此外...

571
来自专栏CDA数据分析师

详解R语言中的遗传算法

前言 人类总是在生活中摸索规律,把规律总结为经验,再把经验传给后人,让后人发现更多的规规律,每一次知识的传递都是一次进化的过程,最终会形成了人类的智慧。自然界规...

23810
来自专栏机器之心

一文读懂遗传算法工作原理(附Python实现)

2895
来自专栏腾讯移动品质中心TMQ的专栏

遗传算法在测试中的应用初探

导读 alphago和master在围棋领域的成绩掀起一股人工智能的热潮之后,人工智能在各个领域的应用成为了大家讨论的焦点。其实机器学习在测试领域的应用也已经有...

2895
来自专栏生信技能树

PGSEA和GSVA你会怎么选择呢?

it tests for each sample whether the average expression of genes in a gene sets ...

1187
来自专栏生信技能树

各种NGS组学数据分析异同点视频讲解

全外显子(Whole-exome sequencing)测序是啥?转录组(RNA-seq)测序是啥?ChIP-seq又是啥?它们之间有什么差别么?傻傻分不清,不...

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

【学习】R语言中的遗传算法

前言 人类总是在生活中摸索规律,把规律总结为经验,再把经验传给后人,让后人发现更多的规规律,每一次知识的传递都是一次进化的过程,最终会形成了人类的智慧。自然界规...

2926
来自专栏锦小年的博客

Nilearn学习笔记1-Nilearn库介绍

nilearn是一个将机器学习、模式识别、多变量分析等技术应用于神经影像数据的应用中,能完成多体素模式分析(MVPA:mutli-voxel pattern a...

2095
来自专栏小鹏的专栏

是AI就躲个飞机-纯Python实现人工智能

代码下载:Here。 很久以前微信流行过一个小游戏:打飞机,这个游戏简单又无聊。在2017年来临之际,我就实现一个超级弱智的人工智能(AI),这货可以躲避从...

8385

扫码关注云+社区