「群体遗传学实战」第三课: 如何对SNP位点进行过滤

往期教程

• 「群体遗传学实战」第一课: 对SNP位点进行注释
• 「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图

---

SNP过滤有两种情况，一种是仅根据位点质量信息（测序深度，回帖质量等）对SNP进行过滤。如果使用GATK对重测序结果进行SNP calling，那么可以考虑下面的标准

• QD< 2.0 || FS> 60.0 || MQ< 40.0 || MQRankSum <−12.5 || ReadPosRankSum <−8.0
• QUAL<30.0||QD<2.0||FS>60.0||MQ<40.0||SOR>4.0和--clusterWindowSize 5 --clusterSize 2

关于这部分的过滤方法，参考如下几篇

• call variant中关于snp筛选的一些思考
• 「简化基因组」如何过滤用GATK分析得到的SNP

另一种过滤会考虑除了测序质量以外的信息，例如文章在方法部分所写的内容

Bi-allelic SNPs with a missing data rate less than 15% and a minor allele count greater than three were kept for population genomic analyses. Additionally, only SNPs at fourfold degenerated sites (89,914 SNPs) were used to construct a neighbor-joining phylogenetic tree using MEGA7 with 500 bootstraps61. ... STRUCTURE analyses were run 20 times for each K value ranging from 2 to 20, using 8,000 randomly selected SNPs at fourfold degenerated sites ...

• Bi-allelic, 相对于multi-allelic, 也就是该位点中只有一个等位基因位点。会过滤掉REF=A, ALT=C,G的SNP位点
• 缺失率低于15%
• 次要等位基因的count数大于3
• 四倍兼并位点

思考题，为什么要用这些规则？

前三个条件的实现相对简单，虽然VCFtools和BCFtools都可以实现这种过滤，但是BCFtools的执行速度更快（大概是前者的2倍），所以我推荐使用BCFtools。

# BCFtools
bcftools view -i 'F_MISSING < 15 & MAC > 3'  -m2 -M2 watermelon_414acc_SNP2.vcf.gz -Oz -o watermelon_414acc_SNP2_flt1.vcf.gz &
# VCFtools
# vcftools --gzvcf watermelon_414acc_SNP2.vcf.gz --min-alleles 2 --max-alleles 2 --max-missing 0.15 --mac 3 --recode --recode-INFO-all  --stdout | bcftools view -Oz -o watermelon_414acc_SNP2_flt1.vcf.gz &
bcftools index watermelon_414acc_SNP2_flt1.vcf.gz

我同时运行了两个程序，最终原始的19,725,853 SNP经BCFtools过滤后为11,925,733，而VCFtools过滤后是12,555,059，BCFtools用时6202秒, VCFtools用时10883秒。我使用vcftools的比较功能，发现问题问题出在MAC的这个标准上，vcftools中--mac 3会包括MAF=3的情况，而我写的bcftools过滤表达式为MAC > 3没有包括3。根据文章的描述，vcftools过滤参数应该写成--mac 4。

出处: Include only sites with Minor Allele Count greater than or equal to the "--mac" value and less than or equal to the "--max-mac" value。

四倍兼并位点(4dTv)过滤稍微麻烦一些，似乎也不是所有文章都会使用该方法。我个人为使用该方法的主要目的是进一步减少SNP的数目，降低后续构建系统发育树和群体结构分析的计算量。

过滤4dTv位点有两种方法，一种是基于注释的VCF文件自己写脚本处理，一种是先生成所有的4DTV候选位置，然后遍历VCF文件并判断当前位点是否为4DTV。此处，我们采用第二种方法，第一种作为练习题。

我们使用Reseqtools根据Fasta和GFF提取所有的4DTV位点

# 提取位点
iTools Fatools getCdsPep -Ref watermelon/97103_genome_v2.fa -Gff watermelon/97103_gene_gff_v2 -4DSite -OutPut watermelon
zcat watermelon.4Dsite.gz | cut -f 1,2 > watermelon.4Dsite.txt

然后我们可以使用BCFtools的-R参数进行过滤，但是速度会很慢，因为每个位点都要和将近400w个位点进行比较。

# 过滤位点
bcftools view -R watermelon.4Dsite.txt watermelon_414acc_SNP2.flt1.vcf.gz -Oz -o watermelon_414acc_SNP2.flt2.vcf.gz

或者我们可以写一个Python脚本，先将所有位置保存在一个集合(set)中，接着遍历VCF文件，将每个位置和存放位置的集合进行比较

python filter_vcf_by_4dtv.py watermelon_414acc_SNP2_flt1.vcf.gz watermelon_414acc_SNP2_flt3.vcf.gz watermelon.4Dsite.txt &

我的脚本运行时间大约是1502s(25分钟)，而用bcftools跑了6小时都还没有结束。

最终19,725,853个SNP经过上述条件过滤后，只剩下了141,324个SNP，和原文的89,914相比，多了大约5万个位点，个人认为是4DTV过滤这一步存在差异。我们之后会用过滤后的位点进行系统发育树构建和群体结构分析。

filter_vcf_by_4dtv.py代码如下