WGS全基因组分析 || SNP过滤

当全基因组分析得到成百上千样本的变异位点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的数据进行过滤分析，一些过滤的手段还比较少遇到。