123456789101112131415161718192021222324252627282930313233 |
- #!/usr/bin/perl
- use strict;
- use warnings;
- my $genome=shift;
- while(<STDIN>){
- 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 '<DUP:TANDEM>' or $alt eq '<DEL>'){
- 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 ;
- }
- }
|