12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758 |
- #!/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 );
- }
|