SNP过滤

SNP过滤

---

SNP过滤

所属目录：紫菜

创建时间：2024/7/25

作者：星云<XingYun>

更新时间：2024/7/25

URL：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));
}

---

总结

本文主要记录了对之前得到的原始的vcf文件进行SNP过滤后，保留过滤后的位点数据。
先用perl脚本统计位点信息，再用R脚本获得过滤后的位点位置信息，最后用vcftools保留过滤后的位点。

–2024/7/25