bcftools/bcf-fix.pl

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4 use warnings;
   5 use Carp;
   6
   7 my $opts = parse_params();
   8 bcf_fix();
   9
  10 exit;
  11
  12 #--------------------------------
  13
  14 sub error
  15 {
  16     my (@msg) = @_;
  17     if ( scalar @msg ) { confess @msg; }
  18     die
  19         "Usage: bcftools view test.bcf | bcf-fix.pl > test.vcf\n",
  20         "Options:\n",
  21         "   -h, -?, --help                  This help message.\n",
  22         "\n";
  23 }
  24
  25
  26 sub parse_params
  27 {
  28     my $opts = {};
  29     while (my $arg=shift(@ARGV))
  30     {
  31         if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }
  32         error("Unknown parameter \"$arg\". Run -h for help.\n");
  33     }
  34     return $opts;
  35 }
  36
  37 sub bcf_fix
  38 {
  39     while (my $line=<STDIN>)
  40     {
  41         if ( $line=~/^#CHROM/ )
  42         {
  43             print
  44 qq[##INFO=<ID=DP4,Number=4,Type=Integer,Description="Read depth for 1) forward reference bases; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref">
  45 ##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for 1) strand bias (exact test); 2) baseQ bias (t-test); 3) mapQ bias (t); 4) tail distance bias (t)">
  46 ##INFO=<ID=AF1,Number=1,Type=Float,Description="EM estimate of site allele frequency without prior">
  47 ##INFO=<ID=AFE,Number=1,Type=Float,Description="Posterior expectation of site allele frequency (with prior)">
  48 ##INFO=<ID=HWE,Number=1,Type=Float,Description="P-value for Hardy-Weinberg equillibrium (chi-square test)">
  49 ##INFO=<ID=MQ,Number=1,Type=Integer,Descriptin="RMS mapping quality">
  50 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
  51 ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
  52 ##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">
  53 ];
  54             print $line;
  55         }
  56         elsif ( $line=~/^#/ )
  57         {
  58             print $line;
  59         }
  60         else
  61         {
  62             my @items = split(/\t/,$line);
  63             my @tags = split(/:/,$items[8]);    # FORMAT tags
  64
  65             my $nidx=2;
  66             my @idxs;   # Mapping which defines new ordering: $idxs[$inew]=$iold; GT comes first, PL second
  67             for (my $i=0; $i<@tags; $i++)
  68             {
  69                 if ( $tags[$i] eq 'GT' ) { $idxs[0]=$i; }
  70                 elsif ( $tags[$i] eq 'PL' ) { $idxs[1]=$i; }
  71                 else { $idxs[$nidx++]=$i; }
  72             }
  73             if ( !exists($tags[0]) or !exists($tags[1]) ) { error("FIXME: expected GT and PL in the format field.\n"); }
  74
  75             # First fix the FORMAT column
  76             $items[8] = 'GT:GL';
  77             for (my $i=2; $i<@tags; $i++)
  78             {
  79                 $items[8] .= ':'.$tags[$idxs[$i]];
  80             }
  81
  82             # Now all the genotype columns
  83             for (my $iitem=9; $iitem<@items; $iitem++)
  84             {
  85                 @tags = split(/:/,$items[$iitem]);
  86                 $items[$iitem] = $tags[$idxs[0]] .':';
  87
  88                 # GL=-PL/10
  89                 my ($a,$b,$c) = split(/,/,$tags[$idxs[1]]);
  90                 $items[$iitem] .= sprintf "%.2f,%.2f,%.2f",-$a/10.,-$b/10.,-$c/10.;
  91
  92                 for (my $itag=2; $itag<@tags; $itag++)
  93                 {
  94                     $items[$iitem] .= ':'.$tags[$idxs[$itag]];
  95                 }
  96             }
  97             print join("\t",@items);
  98         }
  99     }
 100 }
 101