prod_reproduction/202307-IDTWES/script/get_vaf.pl

#!/usr/bin/perl

use strict;
use warnings;
use File::Basename;

my $usage="\n\tperl $0 vcf_file 
\tscript to get chrM/MT VAF infomation.\n\n";
die $usage unless @ARGV;

my $vfile=shift;
my $out=basename($vfile);
$out=~s/.vcf.*//;
my $fh;
if(-T $vfile){
    open $fh,'<',$vfile or die "Cannot open file:$!" ;
}else{
    open $fh, "gzip -dc $vfile |" or die "Cannot open pipe:$! "
}
open my $wfh, '>', $out.'.MT.vaf.xls' or die "Cannot write file: $!\n";

while(<$fh>){
    next if /^#/;
    my ($chr, $pos, $ref, $alt, $format, $call)=(split /\t/)[0,1,3, 4, 8,9];
    # $chr=~s/chr//;
    next unless $chr=~/M/;
    next unless $format=~/AD/;
    next if $call eq '.' or $call eq './.';
    my ($ad_col);
    my @formats=split /:/,$format;
    for my $i (0.. $#formats){
        $ad_col=$i if $formats[$i] eq 'AD';
        # $dp_col=$i if $formats[$i] eq 'DP';
    }
    my ($ad)=(split /:/,$call)[$ad_col];
    if($alt=~/,/){
        my @alts=split /,/, $alt; 
        my @ads=split /,/, $ad;
        my $tot=0;
        for my $cnt(@ads){
            $tot+=$cnt;
        }
        #my $nalt=0;
        next unless $tot>0;
        for my $i(0..$#alts){
            my $base=$alts[$i];
            my $nbase=(split /,/, $ad )[$i+1];
            next if $base=~/\*/;
            print $wfh join("\t", $chr.'-'.$pos, $ref, $base ,100*sprintf("%.2f",$nbase/$tot).'%' ), "\n";
        }
    }else{
        my($nref, $nalt)=split /,/, $ad;              
        print $wfh join("\t", $chr.'-'.$pos, $ref, $alt ,100*sprintf("%.2f",$nalt/($nref+$nalt) ).'%' ), "\n" if $nref+$nalt>0 ;
    }
       
}