misc/sam2vcf.pl

   1 #!/usr/bin/perl -w
   2 #
   3 # VCF specs: http://www.1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcfv3.2
   4
   5 # Contact: pd3@sanger
   6 # Version: 2009-10-08
   7
   8 use strict;
   9 use warnings;
  10 use Carp;
  11
  12 my $opts = parse_params();
  13 do_pileup_to_vcf($opts);
  14
  15 exit;
  16
  17 #---------------
  18
  19 sub error
  20 {
  21     my (@msg) = @_;
  22     if ( scalar @msg ) { croak(@msg); }
  23     die
  24         "Usage: sam2vcf.pl [OPTIONS] < in.pileup > out.vcf\n",
  25         "Options:\n",
  26         "   -r, -refseq <file.fa>            The reference sequence, required when indels are present.\n",
  27         "   -h, -?, --help                   This help message.\n",
  28         "\n";
  29 }
  30
  31
  32 sub parse_params
  33 {
  34     my %opts = ();
  35
  36     $opts{fh_in}  = *STDIN;
  37     $opts{fh_out} = *STDOUT;
  38
  39     while (my $arg=shift(@ARGV))
  40     {
  41         if ( $arg eq '-r' || $arg eq '--refseq' ) { $opts{refseq}=shift(@ARGV); next; }
  42         if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }
  43
  44         error("Unknown parameter \"$arg\". Run -h for help.\n");
  45     }
  46     return \%opts;
  47 }
  48
  49 sub iupac_to_gtype
  50 {
  51     my ($ref,$base) = @_;
  52     my %iupac = (
  53             'K' => ['G','T'],
  54             'M' => ['A','C'],
  55             'S' => ['C','G'],
  56             'R' => ['A','G'],
  57             'W' => ['A','T'],
  58             'Y' => ['C','T'],
  59             );
  60     if ( !exists($iupac{$base}) )
  61     {
  62         if ( $ref eq $base ) { return ('.','0|0'); }
  63         return ($base,'1|1');
  64     }
  65     my $gt = $iupac{$base};
  66     if ( $$gt[0] eq $ref  ) { return ($$gt[1],'0|1'); }
  67     elsif ( $$gt[1] eq $ref ) { return ($$gt[0],'0|1'); }
  68     return ("$$gt[0],$$gt[1]",'1|2');
  69 }
  70
  71
  72 sub parse_indel
  73 {
  74     my ($cons) = @_;
  75     if ( $cons=~/^-/ )
  76     {
  77         my $len = length($');
  78         return "D$len";
  79     }
  80     elsif ( $cons=~/^\+/ ) { return "I$'"; }
  81     elsif ( $cons eq '*' ) { return undef; }
  82     error("FIXME: could not parse [$cons]\n");
  83 }
  84
  85
  86 # An example of the pileup format:
  87 #   1       3000011 C       C       32      0       98      1       ^~,     A
  88 #   1       3002155 *       +T/+T   53      119     52      5       +T      *       4       1       0
  89 #   1       3003094 *       -TT/-TT 31      164     60      11      -TT     *       5       6       0
  90 #   1       3073986 *       */-AAAAAAAAAAAAAA       3       3       45      9       *       -AAAAAAAAAAAAAA 7       2       0
  91 #
  92 sub do_pileup_to_vcf
  93 {
  94     my ($opts) = @_;
  95
  96     my $fh_in  = $$opts{fh_in};
  97     my $fh_out = $$opts{fh_out};
  98     my ($prev_chr,$prev_pos,$prev_ref);
  99     my $refseq;
 100
 101     while (my $line=<$fh_in>)
 102     {
 103         chomp($line);
 104         my ($chr,$pos,$ref,$cons,$cons_qual,$snp_qual,$rms_qual,$depth,@items) = split(/\t/,$line);
 105
 106         my ($alt,$gt);
 107         if ( $ref eq '*' )
 108         {
 109             # An indel is involved.
 110             if ($chr ne $prev_chr || $pos ne $prev_pos)
 111             {
 112                 if ( !$$opts{refseq} ) { error("Cannot do indels without the reference.\n"); }
 113                 if ( !$refseq ) { $refseq = Fasta->new(file=>$$opts{refseq}); }
 114                 $ref = $refseq->get_base($chr,$pos);
 115             }
 116             else { $ref = $prev_ref; }
 117
 118             # One of the alleles can be a reference and it can come in arbitrary order
 119             my ($al1,$al2) = split(m{/},$cons);
 120             my $alt1 = parse_indel($al1);
 121             my $alt2 = parse_indel($al2);
 122             if ( !$alt1 && !$alt2 ) { error("FIXME: could not parse indel:\n", $line); }
 123             if ( $alt1 && $alt2 && $alt1 eq $alt2 ) { $alt2=''; }
 124             if ( !$alt1 )
 125             {
 126                 $alt=$alt2;
 127                 $gt='0|1';
 128             }
 129             elsif ( !$alt2 )
 130             {
 131                 $alt=$alt1;
 132                 $gt='0|1';
 133             }
 134             else
 135             {
 136                 $alt="$alt1,$alt2";
 137                 $gt='1|2';
 138             }
 139         }
 140         else
 141         {
 142             # SNP
 143             ($alt,$gt) = iupac_to_gtype($ref,$cons);
 144         }
 145
 146         print $fh_out "$chr\t$pos\t.\t$ref\t$alt\t$snp_qual\t0\t\tGT:GQ:DP\t$gt:$cons_qual:$depth\n";
 147
 148         $prev_ref = $ref;
 149         $prev_pos = $pos;
 150         $prev_chr = $chr;
 151     }
 152 }
 153
 154
 155 #------------- Fasta --------------------
 156 #
 157 # Uses samtools to get a requested base from a fasta file. For efficiency, preloads
 158 #   a chunk to memory. The size of the cached sequence can be controlled by the 'size'
 159 #   parameter.
 160 #
 161 package Fasta;
 162
 163 use strict;
 164 use warnings;
 165 use Carp;
 166
 167 sub Fasta::new
 168 {
 169     my ($class,@args) = @_;
 170     my $self = @args ? {@args} : {};
 171     if ( !$$self{file} ) { $self->throw(qq[Missing the parameter "file"\n]); }
 172     $$self{chr}  = undef;
 173     $$self{from} = undef;
 174     $$self{to}   = undef;
 175     if ( !$$self{size} ) { $$self{size}=10_000_000; }
 176     bless $self, ref($class) || $class;
 177     return $self;
 178 }
 179
 180 sub read_chunk
 181 {
 182     my ($self,$chr,$pos) = @_;
 183     my $to = $pos + $$self{size};
 184     my $cmd = "samtools faidx $$self{file} $chr:$pos-$to";
 185     my @out = `$cmd`;
 186     if ( $? ) { $self->throw("$cmd: $!"); }
 187     my $line = shift(@out);
 188     if ( !($line=~/^>$chr:(\d+)-(\d+)/) ) { $self->throw("Could not parse: $line"); }
 189     $$self{chr}  = $chr;
 190     $$self{from} = $1;
 191     $$self{to}   = $2;
 192     my $chunk = '';
 193     while ($line=shift(@out))
 194     {
 195         chomp($line);
 196         $chunk .= $line;
 197     }
 198     $$self{chunk} = $chunk;
 199     return;
 200 }
 201
 202 sub get_base
 203 {
 204     my ($self,$chr,$pos) = @_;
 205     if ( !$$self{chr} || $chr ne $$self{chr} || $pos<$$self{from} || $pos>$$self{to} )
 206     {
 207         $self->read_chunk($chr,$pos);
 208     }
 209     my $idx = $pos - $$self{from};
 210     return substr($$self{chunk},$idx,1);
 211 }
 212
 213 sub throw
 214 {
 215     my ($self,@msg) = @_;
 216     croak(@msg);
 217 }