原文地址: SNP Filtering Tutorial
感覺文章寫的真棒，于是就翻譯了一遍
目標(biāo)
- 學(xué)習(xí)如何使用VCFtools根據(jù)缺失信息，基因型深度(genotype depth)，位點(diǎn)質(zhì)量得分，次等位基因頻率和基因型總深度(genotype call depth)去過濾VCF文件
- 學(xué)習(xí)使用vcflib過濾 FreeBayes分析RAD數(shù)據(jù)后得到VCF文件
- 根據(jù)群體內(nèi)HWE過濾VCF
- 學(xué)習(xí)將VCF文件分開輸出為SNP和INDEL
- 使用單倍型鑒定腳本（haplotyping scripts）進(jìn)一步過濾SNP，處理旁系同源基因和基因型識(shí)別時(shí)的報(bào)錯(cuò)
正文
歡迎大家來到SNP過濾練習(xí)部分。我們第一部分的練習(xí)適用于基本上所有的VCF文件。第二部分練習(xí)，我們將會(huì)處理FreeBayes得到的VCF文件。當(dāng)然，其他SNP caller也能輸出類似的注釋信息，所以主要學(xué)習(xí)分析思想，而不是直接套用代碼。
讓我們先來創(chuàng)建分析的工作文件目錄
mkdir filtering
cd filtering
然后大家根據(jù)我提供的微云地址下載后續(xù)需要的數(shù)據(jù)，地址為https://share.weiyun.com/5VeHoQ3 密碼：xqeyqf
unzip data.zip
ls -l
或者你可以嘗試原文提供的下載方式。
通用過濾
在熱身階段，我們會(huì)用VCFtools(http://vcftools.sourceforge.net) 對(duì)我們的vcf文件進(jìn)行一波過濾。該程序由一個(gè)二進(jìn)制運(yùn)行程序和一些perl腳本組成，用途就和名字一樣，就是用于處理VCF。(PS: 作者認(rèn)為0.1.11版本的vcftools有更多使用的過濾命令，而我用的是0.1.17)
我們先處理 raw.vcf , 該文件存在很多有可能出錯(cuò)的位點(diǎn)，并且還有很多位點(diǎn)只在一個(gè)樣本出現(xiàn)。我們將會(huì)分成三步進(jìn)行講解
第一步: 我們只保留那些一般群體中出現(xiàn)的位點(diǎn)，并且最低的質(zhì)量分?jǐn)?shù)不低于30，次等位基因深度(minor allele count)不少于3
vcftools --gzvcf raw.vcf.gz --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out raw.g5mac3
在這行代碼中，我們用--gzvcf傳遞了一個(gè)壓縮的vcf文件，--max-missing 0.5表示過濾掉低于缺失率高于50%的位點(diǎn)，--mac 3則是過濾 次等位基因深度低于3。該參數(shù)設(shè)置和具體的群體有關(guān)，這里用的群體至少包含了1個(gè)純合個(gè)體，1個(gè)或3個(gè)雜合個(gè)體。
--recode表示輸出過濾后的VCF文件，--recode-INFO-all表示包含原來文件中所有INFO信息，--out則是輸出文件的前綴。
代碼運(yùn)行信息可能會(huì)是下面這個(gè)樣子
After filtering, kept 40 out of 40 Individuals
Outputting VCF file...
After filtering, kept 78434 out of a possible 147540 Sites
Run Time = 13.00 seconds
兩個(gè)簡(jiǎn)單的過濾標(biāo)準(zhǔn)就干掉了將近50%的數(shù)據(jù)，這樣子后面一些步驟運(yùn)行速度也會(huì)更快了
vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3
現(xiàn)在，我么有了一個(gè)過濾后的VCF文件，即raw.g5mac3.recode.vcf 。下一步就是根據(jù)最低深度進(jìn)行過濾(還可以加上平均深度)
vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3
這一步我們保留了在至少有3個(gè)read的位點(diǎn)。由于GATK/FreeBayes在群體 SNP calling時(shí)會(huì)同時(shí)分析所有的樣本，因此即便位點(diǎn)的覆蓋度低，可信度也很高。
當(dāng)然作者還擔(dān)心你不相信，于是還提供了一個(gè)腳本幫你評(píng)估可能的錯(cuò)誤率
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/ErrorCount.sh
chmod +x ErrorCount.sh
./ErrorCount.sh raw.g5mac3dp3.recode.vcf
輸出信息如下
This script counts the number of potential genotyping errors due to low read depth
It report a low range, based on a 50% binomial probability of observing the second allele in a heterozygote and a high range based on a 25% probability.
Potential genotyping errors from genotypes from only 1 read range from 0 to 0.0
Potential genotyping errors from genotypes from only 2 reads range from 0 to 0.0
Potential genotyping errors from genotypes from only 3 reads range from 15986 to 53714.22
Potential genotyping errors from genotypes from only 4 reads range from 6230 to 31502.04
Potential genotyping errors from genotypes from only 5 reads range from 2493 to 18914
40 number of individuals and 78434 equals 3137360 total genotypes
Total genotypes not counting missing data 2380094
Total potential error rate is between 0.0103815227466 and 0.0437504821238
SCORCHED EARTH SCENARIO
WHAT IF ALL LOW DEPTH HOMOZYGOTE GENOTYPES ARE ERRORS?????
The total SCORCHED EARTH error rate is 0.129149100834.
目前來看，我們的VCF中基因型由于覆蓋度低于5的最大錯(cuò)誤率低于5%，因此沒有啥可以操心。
下一步則是刪除那些測(cè)序結(jié)果不是很好的樣本，也就是缺失值過多的樣本。先統(tǒng)計(jì)下每個(gè)樣本的缺失比例
vcftools --vcf raw.g5mac3dp3.recode.vcf --missing-indv
cat out.imiss
在結(jié)果中你發(fā)現(xiàn)部分個(gè)體的缺失率高達(dá)99.6%，這個(gè)樣本是必然要被刪掉的。不過在此之前，我們可以用柱狀圖來看下分布
mawk '!/IN/' out.imiss | cut -f5 > totalmissing
gnuplot << \EOF
set terminal dumb size 120, 30
set autoscale
unset label
set title "Histogram of % missing data per individual"
set ylabel "Number of Occurrences"
set xlabel "% of missing data"
#set yr [0:100000]
binwidth=0.01
bin(x,width)=width*floor(x/width) + binwidth/2.0
plot 'totalmissing' using (bin($1,binwidth)):(1.0) smooth freq with boxes
pause -1
EOF
絕大部分的個(gè)體的缺失率低于0.5，因此我們用0.5作為我們的閾值進(jìn)行過濾。
那么問題來了，怎么挑選出缺失率大于50%的個(gè)體呢？下面的代碼是其中一種解決方法
mawk '$5 > 0.5' out.imiss | cut -f1 > lowDP.indv
下一步我們用vcftools根據(jù)提供的樣本編號(hào)進(jìn)行過濾
vcftools --vcf raw.g5mac3dp3.recode.vcf --remove lowDP.indv --recode --recode-INFO-all --out raw.g5mac3dplm
作者提供了自動(dòng)化的過濾腳本方便，用于完成剛才的幾步。這個(gè)腳本會(huì)自動(dòng)推斷出比較好的閾值，不過你可以選擇no然后自己輸入一個(gè)閾值。
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/filter_missing_ind.sh
chmod +x filter_missing_ind.sh
./filter_missing_ind.sh raw.g5mac3dp3.recode.vcf DP3g95maf05
在過濾一些樣本后，我們可以根據(jù)最大缺失值比例，平均深度和次等位基因頻率(MAF)進(jìn)行進(jìn)一步過濾
vcftools --vcf raw.g5mac3dplm.recode.vcf --max-missing 0.95 --maf 0.05 --recode --recode-INFO-all --out DP3g95maf05 --min-meanDP 20
這一步過后，我們得到了12,754個(gè)候選位點(diǎn)。
不過上面是對(duì)群體的所有樣本進(jìn)行過濾。假如你的群體來自于多個(gè)區(qū)域，你想對(duì)不同的群體的樣本分別過濾，那么應(yīng)該怎么做呢？VCFtools工具本身不支持這個(gè)功能，但是我們可以將該任務(wù)進(jìn)行拆分成不同步驟，達(dá)到我們想要的結(jié)果。
你得提供一個(gè)包含樣本來源信息的文本，比如說我們解壓縮后得到的popmap
head popmap
BR_002 BR
BR_004 BR
BR_006 BR
BR_009 BR
BR_013 BR
BR_015 BR
BR_016 BR
BR_021 BR
BR_023 BR
BR_024 BR
然后我們按照第二列的來源信息對(duì)文件進(jìn)行拆分
mawk '$2 == "BR"' popmap > 1.keep && mawk '$2 == "WL"' popmap > 2.keep
下一步，用VCFtools分別估計(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
于是我們就會(huì)得到1.lmiss和2.lmiss
我們將兩個(gè)文件進(jìn)行合并，然后根據(jù)最后一列提取出缺失率大于0.1的樣本.
cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci
利用VCFtools進(jìn)行過濾
vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05
作者用兩個(gè)來源的群體作為例子進(jìn)行講解，實(shí)際情況是我們可能會(huì)遇到超過2個(gè)的情況。作者就寫了一個(gè)腳本進(jìn)行處理
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/pop_missing_filter.sh
chmod +x pop_missing_filter.sh
./pop_missing_filter.sh
