#!/usr/bin/perl use strict; use warnings; use File::Basename; my $ref=shift; my $file=shift; my (%dat); open my $vcffh, "gzip -dc $ref |" or die "Cannot open file:$! \n"; while(<$vcffh>){ next if /^#/; my ($chr, $pos, $id,$ref,$alt)=(split /\t/)[0..4]; my $posid=join("\t",$chr, $pos, $ref,$alt); next if $id eq '.'; $dat{$posid}=$id; } close($vcffh); open my $fh, '<', $file or die "Cannot open file: $! \n"; my $sample=basename($file); $sample=~s/.SNPsite.vcf//; print join("\t","#Sample","Chr","Pos","Ref","Alt","ID","Depth","VAF"),"\n"; while(<$fh>){ chomp; next if /^#/; my ($chr,$pos,$ref, $alts,$ft, $calls )=(split /\t/)[0,1,3,4,8,9]; my @alt=split /,/,$alts; for my $i(0..$#alt){ my $call=$alt[$i]; next if $call eq '*'; my $posid=join("\t",$chr, $pos, $ref,$call); my $id=exists $dat{$posid} ? $dat{$posid} : 'NO_ID'; my ($af,$dp)=get_AF($ft,$calls,$i); print join("\t",$sample,$chr,$pos,$ref,$call,$id,$dp,$af),"\n"; } } close($fh); sub get_AF{ my ($ft,$call,$i)=@_; return (' -',0) unless $call=~/:/; my @formats=split /:/,$ft; my $ad_index; for my $j(0..$#formats){ $ad_index=$j if $formats[$j] eq 'AD'; } my $ads=(split /:/,$call)[$ad_index]; return (' -',0) unless $ads =~/,/; my $tot; for my $num(split /,/,$ads){ $tot+=$num; } my $alt_num=(split /,/,$ads)[$i+1]; return(sprintf("%.2f",$alt_num/$tot),$tot ); }