#!/usr/bin/perl use strict; use warnings; my $genome=shift; while(){ if(/^#/){ print ; next; } my ($chr,$pos,$id,$ref,$alt,$qual, $filter,$other)=split /\t/, $_, 8; next unless $filter eq 'PASS'; next if $other=~/IMPRECISE/; if($alt eq '' or $alt eq ''){ my ($end)=$other=~/END=(\d+)/; my $new_end; if ($end - $pos > 10000){ $new_end+=$pos+10000 }elsif($end < $pos){ next }else{ $new_end=$end } my $seq=`samtools faidx $genome $chr:$pos-$new_end |awk 'NR>1'`; $seq=~s/\s//g; print join("\t",$chr,$pos,$id,$seq,$ref,$qual, $filter,$other); }else{ print ; } }