本文轉載於https://www.jianshu.com/p/e6d5dd774c6e
SNP位點過濾
SNP過濾有兩種情況，一種是僅根據位點質量信息（測序深度，回帖質量等）對SNP進行粗過濾。如果使用GATK對重測序結果進行SNP calling，那么可以考慮下面的標准
QD< 2.0 || FS> 60.0 || MQ< 40.0 || MQRankSum <−12.5 || ReadPosRankSum <−8.0QUAL<30.0||QD<2.0||FS>60.0||MQ<40.0||SOR>4.0和--clusterWindowSize 5 --clusterSize 2
另一種過濾會考慮除了測序質量以外的信息，例如文章在方法部分所寫的內容
Bi-allelic SNPs with a missing data rate less than 15% and a minor allele count greater than three were kept for population genomic analyses. Additionally, only SNPs at fourfold degenerated sites (89,914 SNPs) were used to construct a neighbor-joining phylogenetic tree using MEGA7 with 500 bootstraps61. ... STRUCTURE analyses were run 20 times for each K value ranging from 2 to 20, using 8,000 randomly selected SNPs at fourfold degenerated sites ...
- Bi-allelic, 相對於multi-allelic, 也就是該位點中只有一個等位基因位點。會過濾掉REF=A, ALT=C,G的SNP位點
- 缺失率低於15%: 保證對於任意一個SNP，群體里至少有85%樣本有基因型
- 次要等位基因的count數大於3: 在本文語境中相當於MAF=0.01
- 四倍兼並位點: 在2013年的黃瓜NG中，選擇4D位點原因是它們的選擇壓小，更能反應群體結構(population structure and demography)
前三個條件的實現相對簡單，雖然VCFtools和BCFtools都可以實現這種過濾，但是BCFtools的執行速度更快（大概是前者的2倍），所以我推薦使用BCFtools。
1 # BCFtools 2 bcftools view -i 'F_MISSING < 15 & MAC > 3' -m2 -M2 watermelon_414acc_SNP2.vcf.gz -Oz -o watermelon_414acc_SNP2_flt1.vcf.gz & 3 # VCFtools 4 # vcftools --gzvcf watermelon_414acc_SNP2.vcf.gz --min-alleles 2 --max-alleles 2 --max-missing 0.15 --mac 3 --recode --recode-INFO-all --stdout | bcftools view -Oz -o watermelon_414acc_SNP2_flt1.vcf.gz & 5 bcftools index watermelon_414acc_SNP2_flt1.vcf.gz
我同時運行了兩個程序，最終原始的19,725,853 SNP經BCFtools過濾后為11,925,733，而VCFtools過濾后是12,555,059，BCFtools用時6202秒, VCFtools用時10883秒。我使用
vcftools的比較功能，發現問題問題出在MAC的這個標准上，vcftools中
--mac 3會包括MAF=3的情況，而我寫的bcftools過濾表達式為
MAC > 3沒有包括3。根據文章的描述，vcftools過濾參數應該寫成
--mac 4。
四倍兼並位點(4d)過濾稍微麻煩一些，似乎也不是所有文章都會使用該方法。我個人為使用該方法的主要目的是進一步減少SNP的數目，降低后續構建系統發育樹和群體結構分析的計算量。
過濾4d位點有兩種方法，一種是基於注釋的VCF文件自己寫腳本處理，一種是先生成所有的4D候選位置，然后遍歷VCF文件並判斷當前位點是否為4D。此處，我們采用第二種方法，第一種作為練習題。
我們使用
Reseqtools根據Fasta和GFF提取所有的4D位點
1 # 提取位點 2 iTools Fatools getCdsPep -Ref watermelon/97103_genome_v2.fa -Gff watermelon/97103_gene_gff_v2 -4DSite -OutPut watermelon 3 zcat watermelon.4Dsite.gz | cut -f 1,2 > watermelon.4Dsite.txt
然后我們可以使用BCFtools的-R參數進行過濾，但是速度會很慢，因為每個位點都要和將近400w個位點進行比較。、
1 bcftools view -R watermelon.4Dsite.txt watermelon_414acc_SNP2.flt1.vcf.gz -Oz -o watermelon_414acc_SNP2.flt2.vcf.gz
除了4D位點過濾外，更常見的一種過濾方法是基於LD(連鎖不平衡)對SNP進行過濾，我們這里使用Plink進行數據過濾。
Plink的過濾是基於VCF的ID列，而我們這里的數據的ID列標記為缺失，因此我們需要先用bcftools annotate對位點進行簡單注釋。
1 # 需要注釋位點，增加ID列 2 bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' watermelon_414acc_SNP2_flt2.vcf.gz -Oz -o watermelon_414acc_SNP2_flt2_anno.vcf.gz
接着用Plink的--indep-pairwise 窗口大小 步長 R2篩選位點
1 plink --vcf watermelon_414acc_SNP2_flt2_anno.vcf.gz --const-fid --allow-extra-chr --indep-pairwise 50 10 0.2 --out watermelon_414ac 2 c_SNP2_flt3
最后用plink extract根據"prune.in"從原來的vcf文件中提取信息
1 plink --allow-extra-chr --extract watermelon_414acc_SNP2_flt3.prune.in --make-bed --out watermelon_414acc_SNP2_flt3 --recode vcf-iid --vcf watermelon_414acc_SNP2_flt2_anno.vcf.gz