當(dāng)全基因組分析得到成百上千樣本的變異位點(diǎn)vcf文件時(shí)，即可開(kāi)始進(jìn)行下游的一系列分析。而下游分析的首先工作，就是對(duì)得到的raw vcf文件進(jìn)行篩選過(guò)濾?？梢哉f(shuō)過(guò)濾質(zhì)量的好壞，直接決定了下游一系列分析的準(zhǔn)確性與計(jì)算速度。

本文是根據(jù)SNP Filtering Tutorial文檔進(jìn)行的總結(jié)學(xué)習(xí)

數(shù)據(jù)軟件的準(zhǔn)備

教程提供的數(shù)據(jù)是存放在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

過(guò)濾步驟

STEP-1 常見(jiàn)過(guò)濾

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%的個(gè)體有數(shù)據(jù)，就過(guò)濾掉。保留下來(lái)的variants數(shù)據(jù)就是50%以上都有數(shù)據(jù)。
• --mac 3次要等位基因計(jì)數(shù)為3，過(guò)濾小于3的variants
• --recode 表示要輸出一個(gè)新的vcf文件
• --recode-INFO-all：新文件中保留原vcf文件中的所有INFO信息
• -out：輸出的文件名。

STEP-2 覆蓋深度過(guò)濾

vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3 

• --minDP 3此variant上覆蓋的最小reads數(shù)為3，過(guò)濾低于3條reads的variant

STEP-3 去除高缺失率的樣本

### 查看各個(gè)樣本的缺失率情況
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查看每個(gè)樣本的variant缺失情況。輸出文件"out.imiss"，最后一列即表示樣本中的variants缺失率。
• --remove 根據(jù)lowDP.indv列表去除vcf文件中的高缺失率樣本。

STEP-4 filter by population specific call rate

多個(gè)地點(diǎn)的樣本進(jìn)行地點(diǎn)特異的篩選，需準(zhǔn)備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文件中的不同地點(diǎn)列表
csvtk filter2 popmap -t -f '$2 == "BR"' >1.keep
csvtk filter2 popmap -t -f '$2 == "WL"' >2.keep

### 3. 計(jì)算每個(gè)群體中基因座的缺失數(shù)據(jù)
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. 缺失數(shù)據(jù)0.1閾值的基因座列表
cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci

### 5. 根據(jù)列表進(jìn)行過(guò)濾
vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05

• --keep 1.keep --missing-site 估算每個(gè)群體中基因座的缺失數(shù)據(jù)，生成1.lmiss和2.lmiss文件
• --exclude-positions 根據(jù)列表進(jìn)行過(guò)濾。

STEP-5 根據(jù)等位基因平衡進(jìn)行過(guò)濾

根據(jù)vcffilter進(jìn)行過(guò)濾，期望等位基因的頻率接近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 表示過(guò)濾應(yīng)用于每個(gè)位點(diǎn)。

STEP-6 過(guò)濾兩個(gè)鏈都有reads的位點(diǎn)

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 過(guò)濾其它的一些標(biāo)準(zhǔn)

• 因?yàn)镽AD-seq數(shù)據(jù)非隨機(jī)分布，可刪除質(zhì)量值得分低于深度1/4的基因座

vcffilter -f "QUAL / DP > 0.25" DP3g95p5maf05.fil4.vcf > DP3g95p5maf05.fil5.vcf

• 更為復(fù)雜的過(guò)濾標(biāo)準(zhǔn)

### 1.首先生成每個(gè)位點(diǎn)的覆蓋深度
cut -f8 DP3g95p5maf05.fil5.vcf | grep -oe "DP=[0-9]*" | sed -s 's/DP=//g' > DP3g95p5maf05.fil5.DEPTH

### 2.生成質(zhì)量值列表文件
mawk '!/#/' DP3g95p5maf05.fil5.vcf | cut -f1,2,6 > DP3g95p5maf05.fil5.vcf.loci.qual

### 3. 計(jì)算平均覆蓋度
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上沒(méi)有質(zhì)量分?jǐn)?shù)2倍覆蓋度的位點(diǎn)
paste DP3g95p5maf05.fil5.vcf.loci.qual DP3g95p5maf05.fil5.DEPTH | mawk -v x=2084 '$4 > x' | mawk '$3 < 2 * $4' > DP3g95p5maf05.fil5.lowQDloci

### 6. 刪除指定的位點(diǎn)，并使用vcftools重新計(jì)算loci的覆蓋深度
vcftools --vcf DP3g95p5maf05.fil5.vcf --site-depth --exclude-positions DP3g95p5maf05.fil5.lowQDloci --out DP3g95p5maf05.fil5

### 7. 輸出文件并保留深度分?jǐn)?shù)
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以上的位點(diǎn)
vcftools --vcf  DP3g95p5maf05.fil5.vcf --recode-INFO-all --out DP3g95p5maf05.FIL --max-meanDP 102.5 --exclude-positions DP3g95p5maf05.fil5.lowQDloci --recode 

后記，對(duì)此學(xué)習(xí)教程學(xué)習(xí)過(guò)后學(xué)習(xí)了不少新的SNP過(guò)濾方法，但也得注意此教程主要處理的還是針對(duì)RAD-seq的數(shù)據(jù)進(jìn)行過(guò)濾分析，一些過(guò)濾的手段還比較少遇到。