【直播】 我的基因组27-先简单统计一下全基因组变异情况吧

      在这之前我们已经对VCF格式记录的变异文件有了初步了解，那么接下来我们就实战一下，凭自己的理解来对VCF文件中的内容做一个简单的统计。

      我们提到过VCF文件的第八列是比较复杂的，但是它具体信息在VCF的头文件里面都有详细的说明：

DP=123;VDB=0.498987;SGB=-0.692352;RPB=0.997169;MQB=0.943569;MQSB=0.985669;BQB=0.323648;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=50,35,12,9;MQ=41

      一般我只关心的是DP和DP4的值，下面也是基于此值进行统计。第9,10列是个体的基因型信息，而基因型的推断主要来自于DP和DP4的值。

GT:PL 0/1:182,0,255

      我对于整个VCF文件的行数进行了统计，发现我有480万的变异位点。其中INDEL有70多万，剩余的都是SNV（经常会被混着SNP来称呼，我们先不纠结这个细节）

首先第五列是该位点被bcftools软件找成变异的质量值，我们可以先统计一下质量值的分布。

perl -alne '{next if/^#/;$tmp=int($F[5]/10);$h{$tmp}++}END{print "$_\t$h{$_}" foreach sort {$a <=> $b} keys %h}' realign.vcf >vcf.quality.stat

结果如下：

    我把质量值取了10的区间，上面的0和1代表着质量值低于20的，可以看到质量值低于20的并不多，一般来说，软件call到的variation的质量值低于20，我们会舍弃掉，因为这样的variation可信度不高。

    同理，我们可以统计一下测序深度的分布，每个位点的测序深度在DP=后面的数字，脚本如下：

perl -alne '{next if/^#/; /DP=(\d+);/;$tmp=$1;next unless $tmp;$tmp=int($tmp/10);$h{$tmp}++;}END{print "$_\t$h{$_}" foreach sort {$a <=> $b} keys %h}' realign.vcf

    跟上面脚本的区别是这次无法直接取出某一列，需要捕获DP=后面的数字，同样我也取了10的区间，结果如下：

       可以看到这些变异位点的测序深度分布还是正常的，其中那278469万个测序深度低于10X的变异位点，我们认为置信度不高，后面的分析，可以考虑舍弃掉。

     再来统计一下杂合变异位点跟纯合变异位点的分布吧~

    这次我准备把perl和shell脚本夹着用，对INDEL的统计结果如下：

grep INDEL realign.vcf |perl -alne '{@tmp=split/:/,$F[9];print $tmp[0]}' |sort |uniq -c

因为我用bcftools工具得到的vcf文件里面注释了INDEL与否的信息，所以我可以直接搜索。

结果如下：

335114 0/1

330503 1/1

39404 1/2

杂合纯和比例差不多，说明本次call variation的步骤还算合理。

grep -v INDEL realign.vcf |perl -alne '{@tmp=split/:/,$F[9];print $tmp[0]}' |sort |uniq -c

对SNV的统计结果如下：

2485569 0/1

1645957 1/1

2436 1/2

以上便是对变异的简单统计 。

直播我的基因组的交流群有人问到过这个问题，为什么杂合纯合比例差不多说明本次call variation的步骤还算合理。合理的范围是多少呢？

其实很简单，可以统计一下1000 genome计划里面所有人的vcf文件，看看他们的杂合纯合比例是多少，我们再下一讲就讲这个！

文：Jimmy、阿尔的太阳

图文编辑：吃瓜群众