在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獲取純合與雜合基因型的個數(shù)孝扛；
vcftools --gzvcf output.vcf.gz --hardy --out output.hardy.txt
output.hardy.txt文件的第三列為（HOM1/HET/HOM2），即等位1純合基因型的個數(shù)/雜合基因型的個數(shù)/等位2純合基因型的個數(shù)幽崩。運用perl腳本計算每個位點的缺失率苦始、雜合率及MAF，stat_mr_het_maf.pl:
<>;
my $total = 200;#群體個數(shù)
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.設(shè)置過濾條件
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參數(shù)進行慌申。
可以將這些腳本與命令行整合到一個腳本中陌选，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;#群體個數(shù)
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";
調(diào)用執(zhí)行：
perl filt_vcf.pl input.vcf.gz
