用GATK进行二代测序数据 SNP Calling 流程：（四）变异过滤

GATK推荐的最好的过滤方式是用 VQSR功能，它通过机器学习算法来判断SNP的优劣，因此至少需要两个已存在的 SNP 数据集，一个是经过验证的高质量 SNP 数据集作为真集（如 HapMap），还需要一个质量不是特别高，允许存在小部分假阳性的数据集做训练集（如，1000G）。这些数据集在人类研究中很容易找到，但是在植物中比较困难，因此本流程用硬过滤（hard-filtering）的方法进行变异过滤。

提取SNP和INDEL

SNP 和 INDEL 的过滤参数有所不同，因此分开过滤。

#vcf索引
nohup gatk SelectVariants -V raw.vcf.gz -select-type SNP -O snp..vcf.gz &
nohup gatk SelectVariants -V raw.vcf.gz -select-type INDEL -O indel.vcf.gz &

SNP 过滤

查看过滤参数的分布情况

可以先查看一下想用的过滤参数的分布情况

#使用gatk的variantstotable功能提取过滤参数信息
nohup gatk VariantsToTable -V snp.vcf.gz -F CHROM -F POS -F QD -F QUAL -F SOR -F FS -F MQ -F MQRankSum -F ReadPosRankSum -O snps.recode.table &
nohup gatk VariantsToTable -V indel.vcf.gz -F CHROM -F POS -F QD -F QUAL -F SOR -F FS -F MQ -F MQRankSum -F ReadPosRankSum -O indel.recode.table &

使用gatk VariantFiltration 过滤 SNP

版本：The Genome Analysis Toolkit (GATK) v4.1.8.0、Picard Version: 2.22.8

gatk=gatk
# gz格式需要用tabix索引
vcf=snp.vcf.gz
out_prefix=snp.hardfiltereannotated$gatk VariantFiltration \-V $vcf \-filter "QD < 2.0" --filter-name "QD2" \-filter "QUAL < 30.0" --filter-name "QUAL30" \-filter "FS > 60.0" --filter-name "FS60" \-filter "MQ < 40.0" --filter-name "MQ40" \-filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \-filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \-O ${out_prefix}.vcf

过滤参数使用GATK建议参数。SOR参数在很多文章中没有用到，且与FS参数功能重合，此处也不用了
MQRankSum 和ReadPosRankSum的计算需要至少有一个杂合样本，所以有些位点可能会出现没有这两个INFO的情况，如果用复合过滤式会将这些位点过滤掉，因此不用复合式
gatk输出文件格式会自动识别后缀，为方便后期awk处理，不再输出.gz压缩格式

VariantFiltration 功能仅会在 VCF 文件中对每个变异位点添加过滤注释，因此还需要跟过滤注释提取出通过过滤标准的变异。

#提取通过过滤的变异（PASS），对于没有注释 MQRankSum 和ReadPosRankSum的位点也保留下来
nohup awk '/^#/||$7=="PASS"||$7=="MQRankSum-12.5;ReadPosRankSum-8"||$7=="MQRankSum-12.5"||$7=="ReadPosRankSum-8"' snp.hardfiltereannotated.vcf > snp.hardfiltered.vcf &# 剔除非2个等位基因的位点。对于二倍体物种可以仅保留双等位基因位点。
nohup awk '$5 !~ /([[:alpha:]]|*)+,([[:alpha:]]|*)/{print}' snp.hardfiltered.vcf > snp.hardfiltered.biallelic.vcf &
# 用 VCFtools 分别过滤一份完整度80%和50%的vcf
nohup vcftools --vcf snp_hardfiltered_biallelic.vcf --maf 0.05 --max-missing 0.8 --recode --recode-INFO-all --out snp_hardfiltered_biallelic_maf5_ms80 &nohup vcftools --vcf snp_hardfiltered_biallelic.vcf --maf 0.05 --max-missing 0.5 --recode --recode-INFO-all --out snp_hardfiltered_biallelic_maf5_ms50 &

查看个体杂合、HWE、个体和位点深度

#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --het --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --depth --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --site-depth --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_355_biallelic_maf5_ms50.vcf.recode.vcf --hardy --out snp_hardfiltered_biallelic_maf5_ms50 &

INDEL过滤

筛选出<10bp，双等位基因的 INDEL

nohup gatk SelectVariants -V indel.utx.vcf.gz --O indel.biallelic.vcf --max-indel-size 9 --restrict-alleles-to BIALLELIC &

使用gatk VariantFiltration 过滤 indel

## 过滤注释
nohup gatk VariantFiltration -V indel.biallelic.vcf -O indel.biallelic.10bp.hardfilterannotated.vcf -filter "QD < 2.0" --filter-name "QD2" -filter "QUAL < 30.0" --filter-name "QUAL30" -filter "FS > 200.0" --filter-name "FS200" -filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"# 提取通过过滤的变异
nohup awk '/^#/||$7=="PASS"||$7=="ReadPosRankSum-20"' indel.biallelic.10bp.hardfilterannotated.vcf > indel.biallelic.10bp.hardfiltered.vcf &
# 过滤一份完整度 50% 的SNP
nohup vcftools --vcf indel.biallelic.10bp.hardfiltered.vcf --maf 0.05 --max-missing 0.5 --recode --recode-INFO-all --out indel.biallelic.10bp.hardfiltered.maf5ms50 &

用GATK进行二代测序数据 SNP Calling 流程：（四）变异过滤相关推荐

用GATK进行二代测序数据 SNP Calling 流程：（二）bwa比对和HaplotypeCaller 变异检测
1. 创建基因组索引 bwa index genome.fa 2. 查看read group信息,按read group分组, 比对.合并,生成gvcf 由于数据太多,无法存储过多的中间文件,因此写了 ...
【Sentieon】PacBio HiFi三代测序数据SNP/Indel加速分析
Sentieon软件在二代测序中SNP/Indel变异检测流程已非常成熟,并以其检测准确性高和检测速度快而广受业内人士认可.近日,Sentieon推出了DNAscope LongReads分析流程,深 ...
关于文献中二代测序数据下载（NCBI）的问题
关于文献中二代测序数据下载(NCBI)的问题现在二代测序用于生物学研究非常广泛,大部分文章的序列会上传到Sequence Read Archive(SRA)上,这东西也属于NCBI数据库中的吧,我理 ...
DNA甲基化测序数据的分析流程及相关软件总结
目前检测DNA甲基化的方法众多,主要可以分为以下几类(如表1所示): 图片来源(凡时财等,中国科学: 生命科学,2015) <更多精彩,可关注微信公众号:AIPuFuBio,和大型免费综合生物信 ...
二代测序数据统计分析中为什么是负二项分布？
本文转载自"universebiologygirl",已获授权 1. 转录组数据统计推断的难题在RNA-seq中进行两组间的差异分析是最正常不过的了我们在其它实验中同样会遇到类 ...
网络数据包收发流程(四)：协议栈之packet_type
进入函数netif_receive_skb()后,skb正式开始协议栈之旅. 先上图,协议栈大致过程如下所示: 跟OSI七层模型不同,linux根据包结构对网络进行分层. 比如,arp头和ip头都是紧 ...
二代测序linux软件,二代测序数据分析软件包大全
二代测序数据分析软件包大全 Integrated solutions*CLCbio Genomics Workbench-de novoand reference assembly of Sanger ...
三代测序数据纠错的方法、装置和计算机可读存储介质与流程
三代测序数据纠错的方法.装置和计算机可读存储介质与流程文档序号:15616049发布日期:2018-10-09 21:24 导航: X技术> 最新专利>计算;推算;计数设备的制造及其应用 ...
经典：基因组测序数据从头拼接或组装算法的原理
基因组测序数据的拼接/组装 (图片来源:google) 每一个物种的参考基因组序列(reference genome)的产生都要先通过测序的方法,获得基因组的测序读段(reads),然后再进行从头拼接 ...