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

往期教程

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

---

SNP位点过滤

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

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

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

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%: 保证对于任意一个SNP，群体里至少有85%样本有基因型
• 次要等位基因的count数大于3: 在本文语境中相当于MAF=0.01
• 四倍兼并位点: 在2013年的黄瓜NG中，选择4D位点原因是它们的选择压小，更能反应群体结构(population structure and demography)

前三个条件的实现相对简单，虽然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。

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

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

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

# 提取位点
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脚本filter_vcf_by_4d.py，先将所有位置保存在一个集合中，接着遍历VCF文件，将每个位置和存放位置的集合进行比较，使用方法如下

python filter_vcf_by_4d.py watermelon_414acc_SNP2_flt1.vcf.gz watermelon_414acc_SNP2_flt2.vcf.gz watermelon.4Dsite.txt &

我的脚本运行时间大约是1502s(25分钟)，而用bcftools跑了6小时都还没有结束。虽然时间存在区别，但是两者的输出结果完全相同。

除了4D位点过滤外，更常见的一种过滤方法是基于LD(连锁不平衡)对SNP进行过滤，我们这里使用Plink进行数据过滤。

Plink的过滤是基于VCF的ID列，而我们这里的数据的ID列标记为缺失，因此我们需要先用bcftools annotate对位点进行简单注释。

# 需要注释位点，增加ID列
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' watermelon_414acc_SNP2_flt2.vcf.gz -Oz -o watermelon_414acc_SNP2_flt2_anno.vcf.gz

接着用Plink的--indep-pairwise 窗口大小 步长 R2筛选位点

plink --vcf watermelon_414acc_SNP2_flt2_anno.vcf.gz --const-fid --allow-extra-chr --indep-pairwise 50 10 0.2 --out watermelon_414ac
c_SNP2_flt3

最后用plink extract根据"prune.in"从原来的vcf文件中提取信息

plink --allow-extra-chr --extract watermelon_414acc_SNP2_flt3.prune.in --make-bed --out watermelon_414acc_SNP2_flt3 --recode vcf-iid --vcf watermelon_414acc_SNP2_flt2_anno.vcf.gz

最终，19,725,853个SNP经过四个条件过滤后，剩下了141,324个SNP，和原文的89,914相比，多了大约5万个位点，如果继续根据连锁平衡进行过滤，最后会只剩下10,104个SNP。

filter_vcf_by_4d.py脚本如下