在GWAS分析之前，一般需要对基因型文件进行过滤，包括：缺失率、MAF、杂合率、双等位、仅SNP等。
- 运用bcftools进行进行双等位和仅SNP的过滤；
bcftools view -m2 -M2 -v snps input.vcf.gz -Oz -o output.vcf.gz --threads 20
bcftools支持多线程，速度快；
-m2 -M2 指定输出双等位位点；
-v snps 指定仅输出SNP位点。
- 运用vcftools获取纯合与杂合基因型的个数；
vcftools --gzvcf output.vcf.gz --hardy --out output.hardy.txt
output.hardy.txt文件的第三列为（HOM1/HET/HOM2），即等位1纯合基因型的个数/杂合基因型的个数/等位2纯合基因型的个数。运用perl脚本计算每个位点的缺失率、杂合率及MAF，stat_mr_het_maf.pl:
<>;
my $total = 200;#群体个数
print "Chr\tPos\tMR\tHet\tMAF\n";
while (<>) {
my @a = split/\s+/;
my @b = split/\//,$a[2];
print "$a[0]\t$a[1]\t",($total-$b[0]-$b[1]-$b[2])/$total,"\t",$b[1]/($b[0]+$b[1]+$b[2]),"\t",((sort{$a<=>$b}@b[0,2])[0]*2+$b[1])/(2*($b[0]+$b[1]+$b[2])),"\n";
}
perl stat_mr_het_maf.pl output.hardy.txt > output_mr_het_maf.txt
3.设置过滤条件
filt.pl
<>;
while (<>) {
chomp;
my @a = split/\s+/;
if ($a[2] <= 0.05 && $a[3] <= 0.1 && $a[4] >= 0.1) {
print "$a[0]\t$a[1]\n";
}
}
perl filt.pl output_mr_het_maf.txt > keep_snp_sites.txt
利用vcftools将keep_snp_sites.txt文件中位点取出
vcftools --gzvcf output.vcf.gz --positions keep_snp_sites.txt --recode --recode-INFO-all --out final_out
这一步也可以运用bcftools view的-R参数进行。
可以将这些脚本与命令行整合到一个脚本中，filt_vcf.pl：
system"bcftools view -m2 -M2 -v snps $ARGV[0] -Oz -o $ARGV[0].tmp.vcf.gz --threads 20";
system"vcftools --gzvcf $ARGV[0].tmp.vcf.gz --hardy --out $ARGV[0].tmp.hardy.txt";
open IN,"$ARGV[0].tmp.hardy.txt";
open OUT,">$ARGV[0].tmp.keep.txt";
<IN>;
my $total = 200;#群体个数
while (<IN>) {
my @a = split/\s+/;
my @b = split/\//,$a[2];
my $mr = ($total-$b[0]-$b[1]-$b[2])/$total;
my $het = $b[1]/($b[0]+$b[1]+$b[2]);
my $maf = ((sort{$a<=>$b}@b[0,2])[0]*2+$b[1])/(2*($b[0]+$b[1]+$b[2]));
if ($mr <= 0.05 && $het <= 0.1 && $maf >= 0.1) {
print OUT "$a[0]\t$a[1]\n";
}
}
close IN;
close OUT;
system"vcftools --gzvcf $ARGV[0].tmp.vcf.gz --positions $ARGV[0].tmp.keep.txt --recode --recode-INFO-all --out final_filt_$ARGV[0]";
unlink "$ARGV[0].tmp.vcf.gz" "$ARGV[0].tmp.hardy.txt" "$ARGV[0].tmp.keep.txt";
调用执行：
perl filt_vcf.pl input.vcf.gz
