「群体遗传学实战」第三课: 如何对SNP位点进行过滤

往期教程

「群体遗传学实战」第一课: 对SNP位点进行注释
「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图

SNP过滤有两种情况，一种是仅根据位点质量信息（测序深度，回帖质量等）对SNP进行过滤。如果使用GATK对重测序结果进行SNP calling，那么可以考虑下面的标准

QD< 2.0 || FS> 60.0 || MQ< 40.0 || MQRankSum <−12.5 || ReadPosRankSum <−8.0
QUAL<30.0||QD<2.0||FS>60.0||MQ<40.0||SOR>4.0和--clusterWindowSize 5 --clusterSize 2

关于这部分的过滤方法，参考如下几篇

call variant中关于snp筛选的一些思考
「简化基因组」如何过滤用GATK分析得到的SNP

另一种过滤会考虑除了测序质量以外的信息，例如文章在方法部分所写的内容

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%
次要等位基因的count数大于3
四倍兼并位点

思考题，为什么要用这些规则？

前三个条件的实现相对简单，虽然VCFtools和BCFtools都可以实现这种过滤，但是BCFtools的执行速度更快（大概是前者的2倍），所以我推荐使用BCFtools。

# BCFtools
bcftools view -i 'F_MISSING < 15 & MAC > 3'  -m2 -M2 watermelon_414acc_SNP2.vcf.gz -Oz -o watermelon_414acc_SNP2_flt1.vcf.gz &
# VCFtools
# 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 &
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。

出处: Include only sites with Minor Allele Count greater than or equal to the "--mac" value and less than or equal to the "--max-mac" value。

四倍兼并位点(4dTv)过滤稍微麻烦一些，似乎也不是所有文章都会使用该方法。我个人为使用该方法的主要目的是进一步减少SNP的数目，降低后续构建系统发育树和群体结构分析的计算量。

过滤4dTv位点有两种方法，一种是基于注释的VCF文件自己写脚本处理，一种是先生成所有的4DTV候选位置，然后遍历VCF文件并判断当前位点是否为4DTV。此处，我们采用第二种方法，第一种作为练习题。

我们使用Reseqtools根据Fasta和GFF提取所有的4DTV位点

# 提取位点
iTools Fatools getCdsPep -Ref watermelon/97103_genome_v2.fa -Gff watermelon/97103_gene_gff_v2 -4DSite -OutPut watermelon
zcat watermelon.4Dsite.gz | cut -f 1,2 > watermelon.4Dsite.txt

然后我们可以使用BCFtools的-R参数进行过滤，但是速度会很慢，因为每个位点都要和将近400w个位点进行比较。

# 过滤位点
bcftools view -R watermelon.4Dsite.txt watermelon_414acc_SNP2.flt1.vcf.gz -Oz -o watermelon_414acc_SNP2.flt2.vcf.gz

或者我们可以写一个Python脚本，先将所有位置保存在一个集合(set)中，接着遍历VCF文件，将每个位置和存放位置的集合进行比较

python filter_vcf_by_4dtv.py watermelon_414acc_SNP2_flt1.vcf.gz watermelon_414acc_SNP2_flt3.vcf.gz watermelon.4Dsite.txt &

我的脚本运行时间大约是1502s(25分钟)，而用bcftools跑了6小时都还没有结束。

最终19,725,853个SNP经过上述条件过滤后，只剩下了141,324个SNP，和原文的89,914相比，多了大约5万个位点，个人认为是4DTV过滤这一步存在差异。我们之后会用过滤后的位点进行系统发育树构建和群体结构分析。

filter_vcf_by_4dtv.py代码如下

「群体遗传学实战」第三课: 如何对SNP位点进行过滤相关推荐

「群体遗传学实战」第一课: 对SNP位点进行注释
文章我们用于实战的数据集来自于2019年发表于NG的西瓜文章,它提供了GATK过滤后的SNP数据集用于分析,并且附录提供了完整的样本信息.该SNP数据集包括414个西瓜,2000万个SNP信息,数据 ...
「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图
主成分分析(PCA)是一种线性降维方法,能从纷繁复杂的数据中抽离出关键因素,用来区分不同的样本.这里我们不谈PCA背后的数学原理,只谈哪些软件能够处理数据,我找到了以下三款 Plink: https: ...
「技术综述」有三AI不得不看的技术综述
https://www.toutiao.com/i6715153780863664653/ 文/编辑 | 言有三最近遇到了很多新手来交流,网上资料甚多,筛选有时候是个大问题,一般遇到一个新方向,找技 ...
「Vue实战」武装你的前端项目
1. 接口模块处理 1.1 axios二次封装很基础的部分,已封装好的请跳过.这里的封装是依据 JWT import axios from 'axios'import router from '.. ...
react安装_「React实战」三分钟搭建React开发环境
其实16年的时候就已经接触到React,那个时候也只是入门,时隔多年,工作上一直都没有接触到相关的业务,不知不觉,前端的天也开始渐变,看到了很多招聘要求上都是要求会React,三大框架怎么也得熟悉使 ...
springboot redis 刷新时间_「SpringBoot实战」SpringCache + Redis实现数据缓存
关注我的微信公众号:后端技术漫谈不定期推送关于后端开发.爬虫.算法题.数据结构方面的原创技术文章,以及生活中的逸闻趣事. 我目前是一名后端开发工程师.主要关注后端开发,数据安全,网络爬虫,物联网,边 ...
怎么查询redis缓存的数据_阿里开发十年写出这份「Redis简明教程」+「Redis实战」请你查收...
Redis是啥?用Redis官方的话来说就是: Redis is an open source (BSD licensed), in-memory data structure store, used ...
「项目实战」一文读懂思科网络设备IOS系统
今天给大家带来的小知识是一文读懂思科的IOS系统,相信大家都有了解,但是今天呢给大家把完整的流程梳理出来,这样有助于大家记笔记哦! IOS是被用来传送网络服务并启动网络应用的.Cisco路由器的IOS ...
「Nginx实战」中学到的东西用在面试上，面试官都被怼得哑口无言
Nginx到底是什么? Nginx是一个高性能的HTTP和反向代理web服务器,同时也提供了IMAP/POP3/SMTP服务.Nginx是由伊戈尔·赛索耶夫为俄罗斯访问量第二的http://Rambl ...

「群体遗传学实战」第三课: 如何对SNP位点进行过滤

「群体遗传学实战」第三课: 如何对SNP位点进行过滤相关推荐

最新文章

热门文章