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