查看原文
其他

【直播】我的基因组49:Y染色体的SNV不能用常规流程来找?

2017-01-22 Jimmy 生信技能树

在上一次直播中,我们说到了一个不符合我们的认知的问题。就是我的全基因组测序数据里找到的SNV的纯合杂合比例失衡,这着实让我非常纠结。在朋友圈大量求助中,肿瘤所的朋友非常热心的帮我检查了她手头的几百个外显子测序样品,给了我下面这个表格,我简单的截取一部分。从这张表格中可以看到,女性样本X染色体的纯合杂合比例符合我们的认知。


不过,我更好奇女性样本的Y染色体SNV(虽然理论上女性是不可能有Y染色体的)。而且我真正想看的是男性样本的性染色体,在朋友电脑里面只有sort好的bam文件,没有vcf直接统计。所以我就借了朋友的电脑亲自上阵来统计这些指标,把所有她已有的外显子测序文件批量统计一下:



为了节省时间,我就用了bcftools来做SNP-calling,批量统计的代码如下:

  1. ls /media/cancer_path/*bam |while read id

  2. do

  3. file=$(basename $id )

  4. sample=${file%%.*}

  5. echo $sample

  6. samtools mpileup  -r  X -ugf /media/software/bwa/human_g1k_v37.fasta  $id | bcftools call -vmO z -o $sample.chrX.vcf.gz

  7. samtools mpileup  -r  Y -ugf /media/software/bwa/human_g1k_v37.fasta  $id | bcftools call -vmO z -o $sample.chrY.vcf.gz

  8. echo "chrX"

  9. zcat $sample.chrX.vcf.gz  |perl -alne 'next if /^#/;/DP=(\d+);/;print if $1>20'|grep -v "^#" |cut -f 10|cut -d":" -f 1|sort |uniq -c

  10. echo "chrY"

  11. zcat $sample.chrY.vcf.gz  |perl -alne 'next if /^#/;/DP=(\d+);/;print if $1>20'|grep -v "^#" |cut -f 10|cut -d":" -f 1|sort |uniq -c

  12. done

得到的统计表格我稍微进行了整理了(左边是男性,右边是女性):

假设朋友给我提供的性别与样本对应表格是准确无误的!

那么她提供的样本中:女性的X染色体的杂合数量远高于X的纯合。这合情合理,而且女性没有Y染色体,但是X,Y有同源区域,所以女性样本仍然会有Y染色体的SNV,也符合情理,毕竟比例很小嘛。

而她提供的男性样本数据里面出现我现在全基因组数据结果相同的困惑,明明男性只有一条X和一条Y染色体,那么上面的SNV应该是纯合的,但是这里面都是杂合的多于纯合的。跟我面临的情况一模一样!


对此,我提出了几个假设:

1.就是人类的X,Y染色体同源区域太多了,即使是PE150的建库测序策略也无法保证reads正确的匹配到参考基因组应有的位置。

2.参考基因组在这两条染色体本来就是模糊不清。

3.我们常规的SNV calling流程在,X,Y染色体上面,准确率很有限!


既然我已经在大样本里面验证了这个现象,那么可以暂时排除是公司把我的样本弄错了那个假设啦!



接下来,我就需要详细解释我自己提出的3个假设咯!

       同时在这里向朋友圈给我提出各种建议的朋友表示衷心的感谢!

下面是大家的建议列表的部分摘抄:


不是说男性的就一定都是纯合的,只是男女比例不一样。这在之前的gwas中也可以观察到。甚至有可以导致性别完全误判的个例基因组型。

xy是绝大部分是同源的,这个现象正常。再有看下突变比例分布,0/1什么的说明不了太多问题。选uniq的方法是什么?最差的结果是,你的数据是混lane测的,污染了。

可能1: 女性样品污染;

可能2: 搜索gametologys evolution; 

可能3: chr x link gene tends to be duplicate more .


男性中X,Y上出现0/1的情况主要是同源区域导致,这个可以从这些0/1突变所在区域发现,这些突变强烈富集,主要集中在几个同源区域。但是X,Y上1/1的突变就分布均匀很多了




        对了,有朋友反映用我的samtools和bcftools代码报错,我看了一下,只是因为他们的samtools和bcftools没有升级到最新版,所以给大家提醒一下:

  1. ## Download and install samtools

  2. ## http://samtools.sourceforge.net/

  3. ## http://www.htslib.org/doc/samtools.html

  4. cd ~/biosoft

  5. mkdir samtools &&  cd samtools

  6. wget https://github.com/samtools/samtools/releases/download/1.3.1/samtools-1.3.1.tar.bz2

  7. tar xvfj samtools-1.3.1.tar.bz2

  8. cd samtools-1.3.1

  9. ./configure --prefix=/home/jianmingzeng/biosoft/myBin

  10. make

  11. make install

  12. ~/biosoft/myBin/bin/samtools --help

  13. ~/biosoft/myBin/bin/plot-bamstats --help

  14. cd htslib-1.3.1

  15. ./configure --prefix=/home/jianmingzeng/biosoft/myBin

  16. make

  17. make install

  18. ~/biosoft/myBin/bin/tabix

  19. ## Download and install bcftools

  20. ## http://www.htslib.org/download/

  21. ## http://www.htslib.org/doc/bcftools-1.0.html

  22. cd ~/biosoft

  23. mkdir bcftools && cd bcftools

  24. wget https://github.com/samtools/bcftools/releases/download/1.3.1/bcftools-1.3.1.tar.bz2

  25. tar xvfj bcftools-1.3.1.tar.bz2

  26. cd bcftools-1.3.1

  27. make

  28. cp bcftools /home/jianmingzeng/biosoft/myBin

  29. ~/biosoft/myBin/bin/bcftools --helpniq -c

  30. done


文:Jimmy

图文编辑:吃瓜群众



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

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