当全基因组分析得到成百上千样本的变异位点vcf文件时，即可开始进行下游的一系列分析。而下游分析的首先工作，就是对得到的raw vcf文件进行筛选过滤。可以说过滤质量的好坏，直接决定了下游一系列分析的准确性与计算速度。
本文是根据SNP Filtering Tutorial文档进行的总结学习
数据软件的准备
教程提供的数据是存放在dropbox中，需要梯子。
### 原始下载地址 curl -L -o data.zip https://www.dropbox.com/sh/bf9jxviaoq57s5v/AAD2Kv5SPpHlZ7LC7sBz4va8a?dl=1
unzip data.zip
ls -lh
### vcftools mawk 的下载
conda install mawk vcftools -y
conda install vcflib -y
过滤步骤
STEP-1 常见过滤
vcftools --gzvcf raw.vcf.gz --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out raw.g5mac3
-
--gzvcf输入的vcf.gz文件 -
--max-missing 0.5即missing rate率，某variants在所有样本中低于50%的个体有数据，就过滤掉。保留下来的variants数据就是50%以上都有数据。 -
--mac 3次要等位基因计数为3，过滤小于3的variants -
--recode表示要输出一个新的vcf文件 -
--recode-INFO-all：新文件中保留原vcf文件中的所有INFO信息 -
-out：输出的文件名。
STEP-2 覆盖深度过滤
vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3
-
--minDP 3此variant上覆盖的最小reads数为3，过滤低于3条reads的variant
STEP-3 去除高缺失率的样本
### 查看各个样本的缺失率情况
vcftools --vcf raw.g5mac3dp3.recode.vcf --missing-indv
### 利用csvtk输出缺失率>0.5的样本
csvtk filter2 out.imiss -t -f '$5 > 0.5'|cut -f 1 > lowDP.indv
### 去除高缺失率样本
vcftools --vcf raw.g5mac3dp3.recode.vcf --remove lowDP.indv --recode --recode-INFO-all --out raw.g5mac3dplm
-
--missing-indv查看每个样本的variant缺失情况。输出文件"out.imiss"，最后一列即表示样本中的variants缺失率。 -
--remove根据lowDP.indv列表去除vcf文件中的高缺失率样本。
STEP-4 filter by population specific call rate
多个地点的样本进行地点特异的筛选，需准备popmap文件
### intergrity 0.95，maf 0.05
vcftools --vcf raw.g5mac3dplm.recode.vcf --max-missing 0.95 --maf 0.05 --recode --recode-INFO-all --out DP3g95maf05 --min-meanDP 20
### 2 选择popmap文件中的不同地点列表
csvtk filter2 popmap -t -f '$2 == "BR"' >1.keep
csvtk filter2 popmap -t -f '$2 == "WL"' >2.keep
### 3. 计算每个群体中基因座的缺失数据
vcftools --vcf DP3g95maf05.recode.vcf --keep 1.keep --missing-site --out 1
vcftools --vcf DP3g95maf05.recode.vcf --keep 2.keep --missing-site --out 2
### 4. 缺失数据0.1阈值的基因座列表
cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci
### 5. 根据列表进行过滤
vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05
-
--keep 1.keep --missing-site估算每个群体中基因座的缺失数据，生成1.lmiss和2.lmiss文件 -
--exclude-positions根据列表进行过滤。
STEP-5 根据等位基因平衡进行过滤
根据vcffilter进行过滤，期望等位基因的频率接近0.5
vcffilter -s -f "AB > 0.25 & AB < 0.75 | AB < 0.01" DP3g95p5maf05.recode.vcf > DP3g95p5maf05.fil1.vcf
mawk '!/#/' DP3g95p5maf05.recode.vcf | wc -l
mawk '!/#/' DP3g95p5maf05.fil1.vcf | wc -l
-
AB >0.25 & AB <0.75保留等位基因平衡在0.25~0.75之中 -
-s表示过滤应用于每个位点。
STEP-6 过滤两个链都有reads的位点
vcffilter -f "SAF / SAR > 100 & SRF / SRR > 100 | SAR / SAF > 100 & SRR / SRF > 100" -s DP3g95p5maf05.fil1.vcf > DP3g95p5maf05.fil2.vcf
mawk '!/#/' DP3g95p5maf05.fil2.vcf | wc -l
STEP-7 过滤其它的一些标准
- 因为RAD-seq数据非随机分布，可删除质量值得分低于深度1/4的基因座
vcffilter -f "QUAL / DP > 0.25" DP3g95p5maf05.fil4.vcf > DP3g95p5maf05.fil5.vcf
- 更为复杂的过滤标准
### 1.首先生成每个位点的覆盖深度
cut -f8 DP3g95p5maf05.fil5.vcf | grep -oe "DP=[0-9]*" | sed -s 's/DP=//g' > DP3g95p5maf05.fil5.DEPTH
### 2.生成质量值列表文件
mawk '!/#/' DP3g95p5maf05.fil5.vcf | cut -f1,2,6 > DP3g95p5maf05.fil5.vcf.loci.qual
### 3. 计算平均覆盖度
mawk '{ sum += $1; n++ } END { if (n > 0) print sum / n; }' DP3g95p5maf05.fil5.DEPTH
### 4. 平均值+均值3倍值
python -c "print int(1952+3*(1952**0.5))"
### 5. 找到cutoff上没有质量分数2倍覆盖度的位点
paste DP3g95p5maf05.fil5.vcf.loci.qual DP3g95p5maf05.fil5.DEPTH | mawk -v x=2084 '$4 > x' | mawk '$3 < 2 * $4' > DP3g95p5maf05.fil5.lowQDloci
### 6. 删除指定的位点，并使用vcftools重新计算loci的覆盖深度
vcftools --vcf DP3g95p5maf05.fil5.vcf --site-depth --exclude-positions DP3g95p5maf05.fil5.lowQDloci --out DP3g95p5maf05.fil5
### 7. 输出文件并保留深度分数
cut -f3 DP3g95p5maf05.fil5.ldepth > DP3g95p5maf05.fil5.site.depth
mawk '!/D/' DP3g95p5maf05.fil5.site.depth | mawk -v x=31 '{print $1/x}' > meandepthpersite
### 8.高平均深度的基因座是指示旁系同源物或多拷贝基因座。删除平均深度为102.5以上的位点
vcftools --vcf DP3g95p5maf05.fil5.vcf --recode-INFO-all --out DP3g95p5maf05.FIL --max-meanDP 102.5 --exclude-positions DP3g95p5maf05.fil5.lowQDloci --recode
后记，对此学习教程学习过后学习了不少新的SNP过滤方法，但也得注意此教程主要处理的还是针对RAD-seq的数据进行过滤分析，一些过滤的手段还比较少遇到。
网友评论