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