prod_reproduction/202307-IDTWES/script/trim_MTsite_VAF.pl

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