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