原文地址: SNP Filtering Tutorial
FreeBayes输出结果过滤
下面的教程假设你的结果文件来自于FreeBayes的输出。
FreeBayes的输出信息非常丰富，我们可以根据RAD-seq的特点和freebayes提供的信息，进行更加复杂的过滤操作。先让我们看看FreeBayes到底有哪些信息
mawk '/#/' DP3g95maf05.recode.vcf
下面是之后会用于过滤的INFO 标签
- AB: Allele balance at heterozygous sites: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous
- SAF: Number of alternate observations on the forward strand
- SAR: Number of alternate observations on the reverse strand
- SRF: Number of reference observations on the forward strand
- SRR: Number of reference observations on the reverse strand
- MQM: Mean mapping quality of observed alternate alleles
- MQMR: Mean mapping quality of observed reference alleles
- PAIRED: Proportion of observed alternate alleles which are supported by properly paired read fragments
- PAIREDR: Proportion of observed reference alleles which are supported by properly paired read fragments
Allele balance
第一个过滤标准是Allele balance。 所谓Allele balance就是0和1的比值，0表示和基因组序列相同，1表示和基因组序列不同。由于RAD-seq是对基因组默写特定位点进行测序，因此我们认为Allele balance应该是0.5左右。
作者使用了 vcflib 的vcffilter进行过滤
vcffilter -s -f "AB > 0.25 & AB < 0.75 | AB < 0.01" DP3g95p5maf05.recode.vcf > DP3g95p5maf05.fil1.vcf
上面这条语句保留了AB在0.25和0.75之间，或者AB<0.01的位点。对于AB < 0.01，作者认为这是固定变异(fixed variants), -s表示将这条语句作用于所有位点，而不只是allels。
我们比较过滤前后的变异数
mawk '!/#/' DP3g95p5maf05.recode.vcf | wc -l
mawk '!/#/' DP3g95p5maf05.fil1.vcf | wc -l
单链特有变异
下一步是根据变异在两条链上出现的次数来进行过滤。对于GWAS和RNA-seq而言，这个操作非常棒。因为除非你文库片段特别的短或者非常的长，那么SNP应该只会出现在正链或者反链，举个例子
可以通过根据 FreeBayes提供SAF, SAR, SRF, SRR信息计算后进行过滤
vcffilter -f "SAF / SAR > 100 & SRF / SRR > 100 | SAR / SAF > 100 & SRR / SRF > 100" -s DP3g95p5maf05.fil1.vcf > DP3g95p5maf05.fil2.vcf
我们是利用比值的方法进行选择，因此部分无关reads不会将整个位点给移除。换句话说，上面的命令包保留正链SNP覆盖度反链大100倍的座位，或者是反链比正链大100倍。
mawk '!/#/' DP3g95p5maf05.fil2.vcf | wc -l
9491
这一步只移除了小部分的位点，但是这些位点很有可能是旁系同源基因，微生物污染或者PCR导致的问题。你可以在IGV上检查这些位点，不过可以用samtools tview进行查看
samtools tview BR_006-RG.bam reference.fasta -p E28188_L151
正如你所见，这个区域非常的混乱，很有可能这里出现了超过两个的单倍型(haplotypes)
比对质量
下一步是检查等位基因的各自比对质量.
vcffilter -f "MQM / MQMR > 0.9 & MQM / MQMR < 1.05" DP3g95p5maf05.fil2.vcf > DP3g95p5maf05.fil3.vcf
由于RAD-seq座位(loci)和等位基因(alleles)应该全部起始于相同的位置， 所以两者不会相差很大
mawk '!/#/' DP3g95p5maf05.fil3.vcf | wc -l
9229
这一步大概过滤了不到3%的变异位点，我们可以对其中一处进行检查
samtools tview BR_004-RG.bam reference.fasta -p E20_L173
这里出现了多
正确配对
由于从头组装不一定完美，因此有些座位可能只有未配对的reads比对上。这并不是问题，真正的问题是**所有支持参考等位基因(reference allele)都是配对，而支持可选等位基因(alternative allele)都不配对(paired)，而这才是问题所在
vcffilter -f "PAIRED > 0.05 & PAIREDR > 0.05 & PAIREDR / PAIRED < 1.75 & PAIREDR / PAIRED > 0.25 | PAIRED < 0.05 & PAIREDR < 0.05" -s DP3g95p5maf05.fil3.vcf > DP3g95p5maf05.fil4.vcf
mawk '!/#/' DP3g95p5maf05.fil4.vcf | wc -l
9166
虽然我们得到的位点越来越少了，但是我们的信噪比越来越高了，看其中一个被我们过滤的例子
samtools tview BR_006-RG.bam reference.fasta -p E4407_L138
相对深度的位点质量得分
我们下一步过滤将会检查相对于深度(depth)的座位质量分数(locus quality score).
李恒在他的文章http://arxiv.org/pdf/1404.0929.pdf发现了一些有趣现象，也就是质量分数和位点深度之间的关系, 还有一篇博客也对这个问题进行了讨论（见推荐阅读）。结论是，高覆盖度会导致座位质量分数的膨胀。李恒认为对于大于平均深度加2-3倍的平均深度平方根的位点深度，它的质量分数会是真实位点的2倍，低于该值则是假阳性变异。
这篇教程的作者认为，对于RAD-seq而言可能还是有点保守，因为RAD-seq的reads并非随机分布。作者基于这一想法，进行两次过滤
第一次过滤，移除任何质量分数低于1/4深度的位点
vcffilter -f "QUAL / DP > 0.25" DP3g95p5maf05.fil4.vcf > DP3g95p5maf05.fil5.vcf
第二种过滤（稍显复杂）。第一步是先记录每个位点的深度
cut -f8 DP3g95p5maf05.fil5.vcf | grep -oe "DP=[0-9]*" | sed -s 's/DP=//g' > DP3g95p5maf05.fil5.DEPTH
第二步是创建质量得分列表
mawk '!/#/' DP3g95p5maf05.fil5.vcf | cut -f1,2,6 > DP3g95p5maf05.fil5.vcf.loci.qual
第三步，计算平均深度
mawk '{ sum += $1; n++ } END { if (n > 0) print sum / n; }' DP3g95p5maf05.fil5.DEPTH
1952.82
然后计算均值加上3倍均值平方根
python -c "print int(1952+3*(1952**0.5))"
2084
之后将质量和深度合并，寻找超过阈值但是质量得分并没有比深度高出2倍的位点
paste DP3g95p5maf05.fil5.vcf.loci.qual DP3g95p5maf05.fil5.DEPTH | mawk -v x=2084 '$4 > x' | mawk '$3 < 2 * $4' > DP3g95p5maf05.fil5.lowQDloci
最后，用VCFtools移除这些位点，然后重新计算深度
vcftools --vcf DP3g95p5maf05.fil5.vcf --site-depth --exclude-positions DP3g95p5maf05.fil5.lowQDloci --out DP3g95p5maf05.fil5
现在根据VCFtools输出，只保留深度分数列
cut -f3 DP3g95p5maf05.fil5.ldepth > DP3g95p5maf05.fil5.site.depth
然后计算平均深度，即将上述文件除以31
mawk '!/D/' DP3g95p5maf05.fil5.site.depth | mawk -v x=31 '{print $1/x}' > meandepthpersite
那我们就可以绘制柱状图
gnuplot << \EOF
set terminal dumb size 120, 30
set autoscale
set xrange [10:150]
unset label
set title "Histogram of mean depth per site"
set ylabel "Number of Occurrences"
set xlabel "Mean Depth"
binwidth=1
bin(x,width)=width*floor(x/width) + binwidth/2.0
set xtics 5
plot 'meandepthpersite' using (bin($1,binwidth)):(1.0) smooth freq with boxes
pause -1
EOF
那些高均值深度的位点极可能是旁系同源基因也可能是多拷贝位点。无论是什么情况，这里我们都要把它们给删掉。这里是移除所有平均深度大于102.5的位置
vcftools --vcf DP3g95p5maf05.fil5.vcf --recode-INFO-all --out DP3g95p5maf05.FIL --max-meanDP 102.5 --exclude-positions DP3g95p5maf05.fil5.lowQDloci --recode
最终，VCFtools从9164个位点中保留了8417个位点。当然上面这些步骤只是为了讲解而已，作者已经写好了脚本，自动化的实现这一过滤步骤
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/dDocent_filters
chmod +x dDocent_filters
./dDocent_filters
HWE过滤
下一步是利用HWE进行过滤。李恒发现，HWE是另一种过滤的绝好标准。由于群体结构会导致HWE偏移，因此不能直接对所有样本进行应用。我们需要分群体分析HWE。作者的一个PhD学生写了一个Perl脚本能够做该分析
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/filter_hwe_by_pop.pl
chmod +x filter_hwe_by_pop.pl
让我们按照群体对SNP进行过滤吧。 先用vcflib的另一个工具vcfallelicprimitives将结果转成SNP。
vcfallelicprimitives DP3g95p5maf05.FIL.recode.vcf --keep-info --keep-geno > DP3g95p5maf05.prim.vcf
这将会将原先的结果拆分成SNP和INDEL两个部分，同时保留位点的INFO标签和基因型信息。接着用VCFTools移除其中的indels
vcftools --vcf DP3g95p5maf05.prim.vcf --remove-indels --recode --recode-INFO-all --out SNP.DP3g95p5maf05
我们得到了8379个SNP，然后使用HWE过滤
./filter_hwe_by_pop.pl -v SNP.DP3g95p5maf05.recode.vcf -p popmap -o SNP.DP3g95p5maf05.HWE -h 0.01
注: 通常情况下，我不会用如此高的p值-h 0.01, 这仅仅适用于当前这个案例。一般而言，错误的位点的p值会更低，而且会出现很多群体中。
经过上面的步骤，我们已经得到了完全过滤的VCF文件，这些SNP的可信度比较高。恭喜你，完成了SNP过滤这一步.
下面的内容可以不用继续看了，因为是作者开发的一个脚本的使用小教程，该用具称之为rad_haplotyper。这个脚本依赖于一些perl模块，需要进行额外安装。
curl -L -O https://raw.githubusercontent.com/chollenbeck/rad_haplotyper/master/rad_haplotyper.pl
chmod +x rad_haplotyper.pl
这里我们没有太多时间去看完所有的信息，这个工具目前还在开发之中，之后还会增加很多的参数选项。
举个例子:
# rad_haplotyper.pl -v SNP.DP3g95p5maf05.HWE.recode.vcf -x 40 -mp 1 -u 20 -ml 4 -n -r reference.fasta
Removed 0 loci (0 SNPs) with more than 20 SNPs at a locus
Building haplotypes for BR_024
Building haplotypes for BR_028
Building haplotypes for WL_054
Building haplotypes for BR_016
Building haplotypes for BR_009
Building haplotypes for BR_006
Building haplotypes for BR_041
Building haplotypes for BR_040
Building haplotypes for BR_046
Building haplotypes for BR_031
Building haplotypes for BR_025
Building haplotypes for BR_002
Building haplotypes for WL_058
Building haplotypes for WL_057
Building haplotypes for WL_061
Building haplotypes for WL_069
Building haplotypes for WL_070
Building haplotypes for BR_048
Building haplotypes for WL_031
Building haplotypes for WL_056
Building haplotypes for BR_047
Building haplotypes for WL_079
Building haplotypes for WL_080
Building haplotypes for WL_032
Building haplotypes for WL_071
Building haplotypes for WL_081
Building haplotypes for BR_004
Building haplotypes for BR_021
Building haplotypes for BR_015
Building haplotypes for BR_043
Building haplotypes for WL_066
Filtered 26 loci below missing data cutoff
Filtered 66 possible paralogs
Filtered 17 loci with low coverage or genotyping errors
Filtered 0 loci with an excess of haplotypes
这行脚本发现了209个位点需要被删掉。除了在终端输出外，该脚本还会输出一个stats.out文件.可以根据该文件构建出需要被删除的位点
grep FILTERED stats.out | mawk '!/Complex/' | cut -f1 > loci.to.remove
作者又写了一个脚本根据上面的位点进行过滤
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/remove.bad.hap.loci.sh
chmod +x remove.bad.hap.loci.sh
./remove.bad.hap.loci.sh loci.to.remove SNP.DP3g95p5maf05.HWE.recode.vcf
这一步会直接得到最终的VCF文件，称之为 SNP.DP3g95p5maf05.HWE.filtered.vcf
mawk '!/#/' SNP.DP3g95p5maf05.HWE.filtered.vcf | wc -l
最后我们只得到了7,666个SNP！那么会有多少可能的错误呢？
./ErrorCount.sh SNP.DP3g95p5maf05.HWE.filtered.vcf
...
The total SCORCHED EARTH error rate is 0.0330857386294.
