SNP过滤

SNP过滤

文章目录

SNP过滤
- 前言
- 一. 利用Perl脚本get_vcf_stats.pl统计位点信息
- 二. 利用R脚本149toTZC.2allele.filtered.R画图并获得过滤后的位点位置信息
- 三. 用vcftools保留过滤后的位点
- 四、get_vcf_stats.pl 脚本存放处
总结

SNP过滤

所属目录：紫菜创建时间：2024/7/25作者：星云<XingYun>更新时间：2024/7/25URL：https://blog.csdn.net/2301_78630677/article/details/140733551

前言

之前（也就是博客重测序数据处理得到vcf文件）得到了原始的vcf文件，本文就使用该文件进行SNP过滤，以便后续的分析

也就是本文是对博客【重测序数据处理得到vcf文件】的补充

首先对原来的vcf文件进行过滤，得到一个高质量的vcf文件，用于后续的群体遗传结构分析

原始文件：
在这里插入图片描述

推荐阅读：图文详解 VCF 生信格式 (变异信息)

在这里插入图片描述

一. 利用Perl脚本get_vcf_stats.pl统计位点信息

工作路径：之前创建的5.haitanensis文件夹

首先将之前得到的"149toTZC.2allele.vcf.recode.vcf"文件压缩为gz格式

gzip 149toTZC.2allele.vcf.recode.vcf

使用脚本get_vcf_stats.pl

perl get_vcf_stats.pl 149toTZC.2allele.vcf.recode.vcf.gz >149toTZC.2allele.vcf.recode.vcf.stat.txt

查看

less -SN 149toTZC.2allele.vcf.recode.vcf.stat.txt

在这里插入图片描述

二. 利用R脚本149toTZC.2allele.filtered.R画图并获得过滤后的位点位置信息

这段R代码执行了一系列数据读取、处理和可视化操作，主要目标是对一个VCF统计文件进行质量控制过滤，并保存过滤后SNP的位置信息。

library(readr)
library(ggplot2)
library(dplyr)
data<-read.table("149toTZC.2allele.vcf.recode.vcf.stat.txt", header = T)
head(data)
par(mfrow=c(3,4))
plot(density(data$qual))
#plot(density(data$BQR))
plot(density(na.omit(data$QD)))
plot(density(na.omit(data$MQ)))
plot(density(na.omit(data$FS)))
plot(density(na.omit(data$SOR)))
plot(density(na.omit(data$MQR)))
plot(density(na.omit(data$RPR)))
plot(density(na.omit(data$Inb)))
plot(density(na.omit(data$EH)))
plot(density(na.omit(data$miss_rate)))
plot(density(na.omit(data$mean_depth)))
plot(density(na.omit(data$mean_GQ)))
snp=data
head(snp)
keep.snp<-snp$QD>=5 & !is.na(snp$QD) & snp$MQ>=50 & !is.na(snp$MQ) & snp$FS<10 & !is.na(snp$FS) &snp$SOR < 2 & !is.na(snp$SOR) & snp$MQR > -2.5 & snp$MQR<2.5 & !is.na(snp$MQR) &snp$RPR > -3 & snp$RPR<3 & !is.na(snp$RPR) & snp$Inb>-0.5 & !is.na(snp$Inb) & snp$EH<3 & !is.na(snp$EH) &snp$miss_rate <=0.2 & snp$mean_depth>=5 & snp$mean_depth<80 & snp$mean_GQ>=20
snp.filtered <- snp[keep.snp,]
pos.snp<-snp[keep.snp, c('chr','pos')]
head(pos.snp)
write_tsv(pos.snp,"149toTZC.2allele.filtered.pos.txt")

代码解释：

# 1. 加载包
library(readr) 
library(ggplot2)
library(dplyr)#2. 使用read.table函数读取名为254toCja.2SNP.vcf.stat.txt的文件，假设文件有表头，将其存储在data对象中。
data<-read.table("149toTZC.2allele.vcf.recode.vcf.stat.txt", header = T)#显示data数据框的前几行，检查数据是否正确加载。
head(data)
#3. 设置图形布局为3行4列，准备绘制多个子图。
par(mfrow=c(3,4))#4. 接下来的多行plot(density(...))语句
#绘制data中各个变量的密度图，包括qual、QD、MQ、FS、SOR、MQR、RPR、Inb、EH、miss_rate、mean_depth和mean_GQ。注意，na.omit()用于去除NA值后再计算密度。#5. 将data复制给snp，用于后续操作。
snp=data#6. 定义一个逻辑向量keep.snp，用于标识满足一系列质量控制条件的SNP。条件包括但不限于：QD大于等于5、MQ大于等于50、FS小于10、SOR小于2等。
keep.snp<-snp$QD>=5 & !is.na(snp$QD) & ...#7. 使用keep.snp逻辑向量过滤snp数据框，只保留符合条件的行，生成snp.filtered数据框。
snp.filtered <- snp[keep.snp,]#8. 从过滤后的snp数据框中提取chr（染色体）和pos（位置）两列，生成pos.snp数据框，仅包含过滤后的SNP位置信息。
pos.snp<-snp[keep.snp, c('chr','pos')]#显示pos.snp数据框的前几行，确认过滤和提取操作是否成功。
head(pos.snp)#9. 将pos.snp数据框以制表符分隔的文本格式保存到文件254toTZC.2SNP.filtered.pos.txt中，记录过滤后SNP的染色体和位置信息。
write_tsv(pos.snp,"149toTZC.2allele.filtered.pos.txt")

三. 用vcftools保留过滤后的位点

vcftools --gzvcf 149toTZC.2allele.vcf.recode.vcf.gz --positions 149toTZC.2allele.filtered.pos.txt --recode --recode-INFO-all --out 149toTZC.2allele.filtered

过滤missing rate 超过0.1和maf小于0.01的位点

vcftools --vcf 149toTZC.2allele.filtered.recode.vcf --max-missing 0.9 --maf 0.01 --recode --recode-INFO-all --out 149toTZC.2allele.filtered.missing0.1.maf0.01

四、get_vcf_stats.pl 脚本存放处

get_vcf_stats.pl 脚本

#usage: perl get_vcf_stats.pl Trapa_variant.raw.vcf > Trapa_rawvcf_par.txt#!/usr/bin/perl -w
use strict;my $infile=shift;
# open(IN, "$infile") or die "Cannot open $infile: $!";
open(IN, "gzip -c -d $infile |") or die "Cannot open $infile: $!";
print STDOUT "chr\tpos\ttype\tqual\tfilter\tBQR\tFS\tMQ\tMQR\tQD\tRPR\tSOR\tEH\tInb\tmiss_rate\tmean_depth\tmean_GQ\n";	while(<IN>){chomp;next if(/^#/);my @cols=split/\s+/;my ($chr, $pos, $ref, $alt, $qual, $filter, $info, $format)=@cols[0,1,3,4,5,6,7,8];$qual=0 if($qual eq '.');splice(@cols,0,9);my $type='SNP';if($alt eq '.'){$type='nonref';}elsif(length($ref)>1){ ### only deletion here $type='INDEL';}else{my @alt=split/,/,$alt;$type='MSNP' if($#alt>0);foreach (@alt){if(/\*/ or length($_)>1){ # deletion and insert$type='INDEL';}}}my ($BQR, $EH, $FS, $Inb, $MQ, $MQR, $QD, $RPR, $SOR)=get_info($info,  qw(BaseQRankSum ExcessHet FS InbreedingCoeff MQ MQRankSum QD ReadPosRankSum SOR));
# 	@cols=@cols[(0..7,10..17)];my ($miss_rate, $mean_depth, $mean_GQ)=summary_gt($format, \@cols);print STDOUT "$chr\t$pos\t$type\t$qual\t$filter\t$BQR\t$FS\t$MQ\t$MQR\t$QD\t$RPR\t$SOR\t$EH\t$Inb\t$miss_rate\t$mean_depth\t$mean_GQ\n";	
}
close(IN);sub get_info{my ($info, @fields)=@_;my @results=();foreach my $field (@fields){my ($value)=$info=~/$field=([+|\-|\.|e|0-9]+)/;if(defined $value){push @results, $value;}else{push @results, 'NA';}}return(@results);
}sub summary_gt{my ($format, $arr_ref)=@_;my ($miss, $dp_sum, $GQ_sum)=0;my @format=split/:/,$format;my @index=(undef, undef, undef);foreach (0..$#format){if($format[$_] eq 'GT'){$index[0]=$_;}elsif($format[$_] eq 'DP'){$index[1]=$_;}elsif($format[$_] eq 'GQ' or $format[$_] eq 'RGQ' ){$index[2]=$_;}}foreach my $ind (@$arr_ref){my @info=split/:/, $ind;my ($gt, $dp)=@info[@index[0,1]];my $GQ=0;$GQ=$info[$index[2]] if(defined $index[2]);if($gt=~/\./){$miss++;$dp=0;$GQ=0;	}else{$dp=0 unless(defined $dp and $dp ne '.');$GQ=0 unless(defined $GQ and $GQ ne '.');}die "$format $ind, $gt, $dp, $GQ, $index[0], $index[1], $index[2] ," if($dp =~ /\./ or $GQ =~ /\./);$dp_sum+=$dp;$GQ_sum+=$GQ;}my $num_inds=scalar @$arr_ref;my ($dp_mean, $GQ_mean)=(0,0);unless($num_inds == $miss){$dp_mean=sprintf("%.2f", $dp_sum/($num_inds-$miss));$GQ_mean=sprintf("%.2f", $GQ_sum/($num_inds-$miss));}$miss=sprintf("%.2f", $miss/$num_inds);return(($miss,$dp_mean, $GQ_mean));
}