我的bam文件如下:
4.0G Mar 29 06:18 B_marked_fixed.bam
3.8G Mar 29 13:22 D_marked_fixed.bam
4.5G Mar 29 07:26 T_marked_fixed.bam
其中B是正常组织的WES数据,使用varscan找somatic mutation的时候作为normal,然后对另外两个样本(D和T)计算。
从这个bam文件可以看到这个WES测序深度不够高,可能平均就 50X吧,如果是 200X的WES数据的bam应该是有20G左右文件大小。
了解hg19和hg38参考基因组异同
需要知道hg38这个新版参考基因组到底进步在哪里。(自行搜索咯)
首先看somatic mutation个数
统计得到的统计学显著的somatic mutation个数如下:
278 D_varscan.snp.Somatic.hc
222 T_varscan.snp.Somatic.hc
200 d_varscan.snp.Somatic.hc
174 t_varscan.snp.Somatic.hc
如果只看有可能是somatic mutation个数如下:
1426 D_varscan.snp.Somatic
1375 T_varscan.snp.Somatic
1071 d_varscan.snp.Somatic
1001 t_varscan.snp.Somatic
其中大写字母的文件代表是比对到了hg19,小写字母的文件是我比对到hg38后跑varscan得到的。可以看到,如果是比对到hg38参考基因组的,那么找到的变异位点要稍微少一点点,不过我意识到参考基因组的有一些是非染色体的片段,所以我重新看了看染色体个数分布情况。
hg38hg19chr
hg38hg19chr
10181
8161
8122
8142
593
473
7204
8224
675
695
476
10196
567
2137
458
248
289
1129
71510
31410
61011
4511
71012
91012
1513
0113
1414
2714
2615
2415
9716
41516
21617
91317
2518
1318
1619
161819
2720
7620
71021
11421
1222
3222
13X
2228X
417Y
420Y
104215sum
136276sum
左边的是T样本,右边的是D样本,可以看到,换成hg38这个新版人类的参考基因组之后,找到统计学显著的somatic mutation个数显著减少了。
当然了,仅仅是看个数,意义不大,我们需要仔细分析位点。
然后具体到位点
首先可以借用一系列网页工具:
用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响 比如输入 hg19,13,19447703,C,T 但是一般是看protein-coding基因上面的情况
snp-nexus 网页略微有点复杂
或者把位点当做peaks来注释:http://52.32.26.75:3838/peaks_annotation/
可以使用homer来进行注释
其实如果这个位点位于dbSNP数据库,那么接下来一切查询都可以基于rs ID号来进行关联,虽然 rs ID号 也会有些微变化。
因为具体到位点,就涉及到课题组信息了,不便公布,但是思路给大家了,可以是坐标转换,或者以 rs ID号 进行关联比较。最终其实要载入IGV去一对一比较,而且varscan软件给的high confidence的somatic mutation也需要注意,它默认P值卡的是0.05,其实一刀切并不好。
更多
以上我仅仅是比较了在50X这个测序深度下,VARSCAN软件基于不同参考基因组版本的表现问题。
还可以探索不同的软件,或者不同的测序深度。
我这里只是想说,对配对的WES数据来说,找somatic mutation这件事,值得仔细检查,假阳性问题比较严重。
测序深度太低的数据,找somatic突变真是头疼
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请
点击举报。