misc/sam2vcf.pl

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