群体遗传专题：关于SNP的过滤（1）：如何使用vcftools进行SNP过滤

前言

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

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

预备工作

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

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

下载并且安装VCFTOOLS

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

小试牛刀

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

解读一下这个代码做了什么：
在这段代码中，我们调用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中会有这样的信息：

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

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

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

运行完的log file 结果如下：

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

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

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

下载好写好的脚本

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

运行结果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

下周继续更新这个教程的下半部分。有兴趣的朋友想超前学多点的朋友，也可以到原文链接进行学习。

群体遗传专题-重测序技术介绍

群体遗传专题：系统发育基础知识介绍

群体遗传专题-structure图的构建与理解

群体遗传专题-学习运用Fst

群体遗传专题：读取数据和数据的基本格式

群体遗传专题-PCA分析

群体遗传专题：基础知识概念理解

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。