前言

很多关注我这个专题的朋友都在后台留言，说想了解更多关于SNP过滤的东西。是的，确实这是一个比较重要的问题。第一，过滤到一些低质量的SNP可以防止calling一些假阳性的SNP，这些假阳性的SNP会很大程度影响到后续的一系列的分析，例如GWAS等的分析，最后影响相关生物学问题的解答。第二，如果你有很多的个体，往往你的call完SNP后，VCF文件的大小的会比较大，如果不经过过滤，对下游的群体结构分析，PCA等相关的计算，都会对计算机产生巨大的负担（运行速度还有需要的内存等）。讲了那么多，在接下来几周会从简单开始，讨论学习一下SNP的过滤的流程，常用工具，还有通常推荐使用的一些参数。

文章主要是根据网上SNP Filtering Tutorial的课程总结翻译而成，推荐新手跟着一步一步来操作，进而理解掌握。老司机就自行忽略看看就好了。

预备工作

创建文件夹，下载测试数据包，解压文件

mkdir filtering
cd filtering
curl -L -o data.zip https://www.dropbox.com/sh/bf9jxviaoq57s5v/AAD2Kv5SPpHlZ7LC7sBz4va8a?dl=1
unzip data.zip 
ls -thrl 

看看数据包里包含了什么文件

total 73742
-rwxr--r-- 1 hhu pawsey0149   109633 Mar  6  2015 BR_004-RG.bam
-rwxr--r-- 1 hhu pawsey0149   247496 Mar  6  2015 BR_004-RG.bam.bai
-rwxr--r-- 1 hhu pawsey0149   120045 Mar  6  2015 BR_006-RG.bam
-rwxr--r-- 1 hhu pawsey0149   247712 Mar  6  2015 BR_006-RG.bam.bai
drwxr-s--- 2 hhu pawsey0149     4096 Jul  1 14:54 __MACOSX
-rwxr--r-- 1 hhu pawsey0149      399 Mar  6  2015 popmap
-rwxr--r-- 1 hhu pawsey0149 68264393 Mar  6  2015 raw.vcf.gz
-rwxr--r-- 1 hhu pawsey0149  6804314 Mar  6  2015 reference.fasta
-rwxr--r-- 1 hhu pawsey0149  1085387 Mar  6  2015 reference.fasta.amb
-rwxr--r-- 1 hhu pawsey0149  1068884 Mar  6  2015 reference.fasta.ann
-rwxr--r-- 1 hhu pawsey0149  6379720 Mar  6  2015 reference.fasta.bwt
-rwxr--r-- 1 hhu pawsey0149   353544 Mar  6  2015 reference.fasta.clstr
-rwxr--r-- 1 hhu pawsey0149   976388 Mar  6  2015 reference.fasta.fai
-rwxr--r-- 1 hhu pawsey0149  1594913 Mar  6  2015 reference.fasta.pac
-rwxr--r-- 1 hhu pawsey0149  3189872 Mar  6  2015 reference.fasta.sa
-rwxr--r-- 1 hhu pawsey0149   137209 Mar  6  2015 stats.out

下载并且安装VCFTOOLS

git clone https://github.com/vcftools/vcftools.git
cd vcftools
./autogen.sh
./configure
make
make install

该文章也会用到一个叫做mawk的软件，它其实就是更加简单版的awk，也需要下载，我这里偷懒用conda下载，其实上面的vcftools 也可以用conda下载

conda install mawk

小试牛刀

首先做一个过滤，
我们将只保留已那些变异，至少50%的个体成功被snp calling的位点，最低质量分数为30的位点，次要等位基因计数为3的位点。

vcftools --gzvcf raw.vcf.gz --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out raw.g5mac3

解读一下这个代码做了什么：
在这段代码中，我们调用vcftools，在--vcf标志后给它一个vcf文件， - max-missing 0.5告诉它过滤低于50％的基因型（跨所有个体）--mac 3标志告诉它过滤次要等位基因计数小于3的SNP。
这与基因型有关，因此必须在至少1个纯合子和1个杂合子或3个杂合子中的变异调用。 --recode标志告诉程序使用过滤器写入一个新的vcf文件，--recode-INFO-all保留旧vcf文件中的所有INFO标志。--out 输出的文件名称

运行完后，在屏幕或者生成的log file中会有这样的信息：

After filtering, kept 40 out of 40 Individuals
Outputting VCF file...
After filtering, kept 78434 out of a possible 147540 Sites
Run Time = 17.00 seconds

这个output给了你一些基本的信息，这里你所用到的40个个体，全部被保留下来，然后147540中的78434个变异位点经过过滤后被保留下来（几乎一半的变异位点被过滤掉）。运行时间：17秒。

在得到这个过滤完的文件后，继续进一步过滤：

我们接下来过滤的标准是基因型调用的最小深度和最小平均深度。

vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3 

运行完的log file 结果如下：

After filtering, kept 40 out of 40 Individuals
Outputting VCF file...
After filtering, kept 78434 out of a possible 78434 Sites
Run Time = 14.00 seconds

可以惊讶的发现，所有的位点都被保留了下来。这说明，上述第一步得到的过滤的文件中的所有位点，已经满足了最小深度至少大于3的过滤条件。

在这里你可能会提出疑问，最小深度大于3这个条件是否足够严谨filter出置信的SNP （理论上是越高越好，但是过高也不太好，这个内容会再以后的文章再具体讨论）。简短的回答是像FreeBayes和GATK这样复杂的多样本SNP calling的工具可以自信地调用几乎没有读数的基因型，因为所有样本的变异是同时评估预测的。因此，这里的过滤是基于来自所有个体的所有读数的三个reads和还有其他的先验信息。但是，我们还会做很多其他的过滤步骤！

如果你还是担心深度不够的影响，我们也可以测测低深度可能产生的SNP错误率。

下载好写好的脚本

curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/ErrorCount.sh
chmod +x ErrorCount.sh 
./ErrorCount.sh raw.g5mac3dp3.recode.vcf 

该脚本计算由于低读取深度而导致的潜在基因分型错误的数量，它会报告基于观察杂合子中第二等位基因的50％二项概率和基于25％概率的高范围。

运行结果

Potential genotyping errors from genotypes from only 1 read range from 0 to 0.0
Potential genotyping errors from genotypes from only 2 reads range from 0 to 0.0
Potential genotyping errors from genotypes from only 3 reads range from 15986 to 53714.22
Potential genotyping errors from genotypes from only 4 reads range from 6230 to 31502.04
Potential genotyping errors from genotypes from only 5 reads range from 2493 to 18914
40 number of individuals and 78434 equals 3137360 total genotypes
Total genotypes not counting missing data 2380094
Total potential error rate is between 0.0103815227466 and 0.0437504821238
SCORCHED EARTH SCENARIO
WHAT IF ALL LOW DEPTH HOMOZYGOTE GENOTYPES ARE ERRORS?????
The total SCORCHED EARTH error rate is 0.129149100834.

现在，由于基因型少于5的深度，VCF文件的最大错误率小于5％。看，没什么好担心的。

下一步，我们可以去除那些没有被测序的很好的个体，我们可以通过估计个体间数据丢失的水平来过滤。

vcftools --vcf raw.g5mac3dp3.recode.vcf --missing-indv

这会生成一个out.imiss的文件，看看他里面是什么：

cat out.imiss

INDV    N_DATA  N_GENOTYPES_FILTERED    N_MISS  F_MISS
BR_002  78434   0   13063   0.166548
BR_004  78434   0   16084   0.205064
BR_006  78434   0   25029   0.319109
BR_009  78434   0   30481   0.38862
BR_013  78434   0   69317   0.883762
BR_015  78434   0   8861    0.112974
BR_016  78434   0   29789   0.379797
BR_021  78434   0   17422   0.222123
BR_023  78434   0   43913   0.559872
BR_024  78434   0   24220   0.308795
BR_025  78434   0   21998   0.280465
BR_028  78434   0   26786   0.34151
BR_030  78434   0   74724   0.952699
BR_031  78434   0   26488   0.337711
BR_040  78434   0   19492   0.248515
BR_041  78434   0   17107   0.218107
BR_043  78434   0   16384   0.208889
BR_046  78434   0   28770   0.366805
BR_047  78434   0   13258   0.169034
BR_048  78434   0   24505   0.312428
WL_031  78434   0   22566   0.287707
WL_032  78434   0   22604   0.288191
WL_054  78434   0   32902   0.419486
WL_056  78434   0   34106   0.434837
WL_057  78434   0   37556   0.478823
WL_058  78434   0   31448   0.400949
WL_061  78434   0   35671   0.45479
WL_064  78434   0   47816   0.609634
WL_066  78434   0   10062   0.128286
WL_067  78434   0   47940   0.611215
WL_069  78434   0   38260   0.487799
WL_070  78434   0   21188   0.270138
WL_071  78434   0   16692   0.212816
WL_072  78434   0   46347   0.590904
WL_076  78434   0   78178   0.996736
WL_077  78434   0   55193   0.703687
WL_078  78434   0   54400   0.693577
WL_079  78434   0   19457   0.248068
WL_080  78434   0   30076   0.383456
WL_081  78434   0   30334   0.386746

您可以看到有些人的数据丢失率高达99.6％。我们绝对想要过滤掉这些。现在我们需要创建一个缺失数据超过50％的个体的列表。

mawk '$5 > 0.5' out.imiss | cut -f1 > lowDP.indv

现在我们有一个要删除的个体列表，我们可以直接将其提供给VCFtools进行过滤。

vcftools --vcf raw.g5mac3dp3.recode.vcf --remove lowDP.indv --recode --recode-INFO-all --out raw.g5mac3dplm

然后该文章的原作者非常贴心的，

我们可以将变异限制为在高百分比的个体中，并通过基因型的平均深度进行过滤。

vcftools --vcf raw.g5mac3dplm.recode.vcf --max-missing 0.95 --maf 0.05 --recode --recode-INFO-all --out DP3g95maf05 --min-meanDP 20

这对所有个体应用的calling rate 为（95％）。这对群体是采集自两个地方的变异来说是已足够。但是当你有多个采样地点时，你还希望按个体特定的calling rate进行过滤。 VCFtools不会直接计算，但这是一个简单的解决方法。首先，我们需要一个文件来定义地点（个体）。大多数程序希望文件具有两个制表符分隔列。首先是样本名称，第二个是个体的群体类型分配。

cat popmap

BR_002  BR
BR_004  BR
BR_006  BR
BR_009  BR
BR_013  BR
BR_015  BR
BR_016  BR
BR_021  BR
BR_023  BR
BR_024  BR
BR_025  BR
BR_028  BR
BR_030  BR
BR_031  BR
BR_040  BR
BR_041  BR
BR_043  BR
BR_046  BR
BR_047  BR
BR_048  BR
WL_031  WL
WL_032  WL
WL_054  WL
WL_056  WL
WL_057  WL
WL_058  WL
WL_061  WL
WL_064  WL
WL_066  WL
WL_067  WL
WL_069  WL
WL_070  WL
WL_071  WL
WL_072  WL
WL_076  WL
WL_077  WL
WL_078  WL
WL_079  WL
WL_080  WL
WL_081  WL

现在我们需要创建两个列表，这些列表只包含每个个体的个别名称。

mawk '$2 == "BR"' popmap > 1.keep && mawk '$2 == "WL"' popmap > 2.keep

接下来，我们使用VCFtools来估计每个群体中基因座的缺失数据。

vcftools --vcf DP3g95maf05.recode.vcf --keep 1.keep --missing-site --out 1
vcftools --vcf DP3g95maf05.recode.vcf --keep 2.keep --missing-site --out 2 

这会生成两个文件分别是 1.lmiss 和2.lmiss

快速看看文件里面的内容：

head -3 1.lmiss

CHR     POS     N_DATA  N_GENOTYPE_FILTERED     N_MISS  F_MISS
E1_L101 9       34      0                       0       0
E1_L101 15      34      0                       0       0

这里文章的原作者在第一行加了header 让你更容易去理解每一列。但我们感兴趣的是最后一列是该基因座缺失数据的百分比。我们可以合并这两个文件，并制作一个关于要删除的10％缺失数据阈值的基因座列表。请注意，这是丢失数据总体速率的两倍。

cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci

然后，我们将此文件反馈到VCFtools以删除任何基因座。

vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05

同样，我们只有两个个体，所以我们的过滤可以过滤这些。然而，在多地点采样的研究就不能像这样处理了。