Mercurial > repos > saskia-hiltemann > cgatools_v17

comparison tools/cgatools17/testvariants2VCF-v2.pl @ 1:3a2e0f376f26 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Minor change to tv2vcf.xml to allow for workflow automation
author dgdekoning
date Wed, 21 Oct 2015 10:09:15 -0400
parents
children
comparison
equal deleted inserted replaced
0:751b62d30ae1 1:3a2e0f376f26
1 #!/usr/bin/perl -w
2 use strict;
3
4 #Converts a cgatools testvariant output to a multi-sample VCF file (VCF spec 4.1)
5 #Requires cgatools to access reference genome encoded in .crr format (hg18.crr or hg19.crr)
6 #make sure cgatools is in the $PATH
7 #
8 #Example testvariant file:
9 #variantId chromosome begin end varType reference alleleSeq xRef GS19238 GS19239 GS19240
10 #6874944 chr5 20584031 20584032 snp C T dbsnp.119:rs10037487 00 01 11
11 #6874945 chr5 20584031 20584032 sub C TA 01 01 00
12 #6874946 chr5 20584031 20584032 sub C TG 00 01 00
13 #After converting to VCF:
14 #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT GS19238 GS19239 GS19240
15 #chr5 20584031 6874944;6874945;6874946 AC AT,ATA,ATG . PASS NS=51;RSID=dbsnp.119:rs10037487,.,.;AF=0.57,0.02,0.00;DB GT 0/2 ./. 1/1
16 #
17 #Variants that share the same location (chr,begin,end) will be merged into one locus and their flags (0,1,N) will be converted into genotype calls
18 #Samples that are positive for more than two alleles within the same locus will be flagged and their genotype calls set to unknown (./.)
19 #Look at the sample GS19239 as such an example
20 #
21 #For non-SNP locus, VCF requires an extra reference base immediately upstream of the variant locus be included in the REF and ALT columns
22 #
23
24
25 die "Usage: $0 testvarOutput.txt hg19.crr > vcf.txt 2> runlog.txt\n\t" .
26 "testvarOutput.txt = output file from cgatools testvariants\n\t" .
27 "hg19.crr = reference genome in .crr format (must be the same as used in testvariants)\n\t" .
28 "vcf.txt = converted file in vcf format\n\t" .
29 "runlog.txt = log file\n" unless ($#ARGV == 1);
30 die "Fail to find input file " . $ARGV[0] . "\n" unless (-e $ARGV[0]);
31 die "Fail to find .crr file " . $ARGV[1] . "\n" unless (-e $ARGV[1]);
32 die "Fail to execute cgatools: $!\n" if (system ('cgatools > /dev/null 2>&1') != 0);
33
34 my (undef, undef, undef, $mday, $mon, $year) = localtime;
35 my($timestamp)=sprintf ('%s%02d%02d', $year+1900, $mon+1, $mday);
36 print <<EOF;
37 ##fileformat=VCFv4.1
38 ##fileDate=$timestamp
39 ##source=CGA Tools v1.7.1 listvariants/testvariants
40 ##reference=$ARGV[1]
41 ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples with Fully Called Data">
42 ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
43 ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP Membership">
44 ##INFO=<ID=RSID,Number=A,Type=String,Description="dbSNP rs ids">
45 ##FILTER=<ID=s50,Description="Less than 50% of samples are fully called">
46 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
47 EOF
48
49 &testvar2VCF (@ARGV);
50
51 sub testvar2VCF {
52 #variantId=0 chromosome=1 begin=2 end=3 varType=4 reference=5 alleleSeq=6 xRef=7 GS12885-1100-37-ASM=8
53 my ($inFile, $crrFile) = @_;
54 my (@samples, %seen);
55 open (IN, "$inFile") or die "Fail to open input file $inFile\n";
56 while (<IN>) {
57 chomp;
58 my (@fs) = split (/\t/);
59 if (/^variantId/) {
60 @samples = @fs[8..$#fs];
61 print join ("\t", ('#CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', @samples)), "\n";
62 next;
63 }
64 next unless (/^\d/);
65 my ($key) = $fs[1] . '-' . $fs[2] . '-' . $fs[3];
66 if (!%seen || exists $seen{$key}) {
67 push (@{$seen{$key}}, [@fs]);
68 next;
69 }
70 #has a new locus, need to process saved locus
71 &doSavedLocus (\%seen, \@samples, $crrFile);
72 #reinitialize with the new locus
73 %seen = ();
74 push (@{$seen{$key}}, [@fs]);
75 }
76 close IN;
77 &doSavedLocus (\%seen, \@samples, $crrFile); #do the last locus
78 }
79
80 sub doSavedLocus {
81 my ($seen, $samples, $crrFile) = @_;
82 my ($lines) = values %$seen; #get all the lines that belong to the same locus
83 my ($chr, $begin, $end) = ($lines->[0]->[1], $lines->[0]->[2], $lines->[0]->[3]);
84 my ($info_NS) = scalar @$samples; #num. of samples fully called
85 die "no_samples,$chr-$begin-$end\n" unless ($info_NS > 0);
86 my ($info_DB) = ''; #dbSNP memebership
87 my (@info_RSid) = (); #dbSNP rs ids
88 my (%ref, @alt, @id, %alleleCalls);
89 my ($format) = 'GT'; #GT (genotype) is the only value currently populated for each sample
90 my ($qual, $filter, $i) = ('.', 'PASS', 0); #no quality scores provided in the testvar output
91 my (@info_AF);
92 #initialize every allele count to 0
93 for ($i=0; $i<=scalar @$lines; $i++) {$info_AF[$i] = 0;}
94 my ($pos) = $begin;
95 #check the vartypes
96 my (%varTypes) = map {$_->[4] => 1} @$lines; #column 4 in the testvar output is varType
97 my (@varTypeKeys) = keys %varTypes;
98 #check if this locus is SNP only
99 my ($snpOnly) = 0;
100 $snpOnly = 1 if ($#varTypeKeys == 0 && $varTypeKeys[0] eq 'snp');
101 foreach my $line (@$lines) {
102 my ($refBase) = $line->[5];
103 my ($altBase) = $line->[6];
104 if (!$snpOnly) { #not a SNP only locus
105 my ($decodecrrCmd) = join ('', ("cgatools decodecrr --reference $crrFile --range $chr:", $begin - 1 , "-$begin 2> /dev/null"));
106 my ($re) = `$decodecrrCmd`;
107 die "cgatools decodecrr failed $decodecrrCmd ", join ("\t", @$line), "\n" unless (defined $re);
108 chomp ($re);
109 $refBase = $re . $refBase;
110 $altBase = $re . $altBase;
111 $pos = $begin;
112 }
113 $ref{$refBase} = 1;
114 push (@alt, $altBase);
115 push (@id, $line->[0]);
116 if ($line->[7] =~ /dbsnp/) {
117 $info_DB = 'DB';
118 $line->[7] =~ s/;/-/g; #testvar use ; as delimiter which is reserved by the INFO column
119 push (@info_RSid, $line->[7]);
120 } else {
121 push (@info_RSid, '.');
122 }
123 #collect testvar's allele flags (00, 01, 11, 0N, NN, ...), which starts at column 8.
124 for ($i=0; $i<scalar @$samples; $i++) {
125 push (@{$alleleCalls{$samples->[$i]}}, $line->[$i+8]);
126 }
127 }
128 $pos = $begin + 1 if ($snpOnly);
129 my (@refBases) = keys %ref;
130 die "mismatched_ref_base, " . join (",", (@id, @refBases)) . "\n" if ($#refBases != 0);
131 my ($idstr) = join (";", @id);
132 my ($altstr) = join (",", @alt);
133 print join ("\t", ($chr, $pos, $idstr, $refBases[0], $altstr,$qual)), "\t";
134 my (@allGTcalls) = ();
135 #convert allele flags to genotype calls
136 foreach my $sample (@$samples) {
137 my ($sampleAlleleCalls) = $alleleCalls{$sample};
138 my ($sampleGTcall);
139 if ($sampleAlleleCalls->[0] =~ /\S\S/) { #diploid site
140 $sampleGTcall = &mkDiploidGTcall ($sampleAlleleCalls);
141 } else {
142 $sampleGTcall = &mkHaploidGTcall ($sampleAlleleCalls);
143 }
144 if ($sampleGTcall =~ /warning/) {
145 print STDERR "$sampleGTcall,$chr-$begin-$end,$idstr,$altstr,$sample,", join ("-", @$sampleAlleleCalls), "\n";
146 $sampleGTcall = '.';
147 $sampleGTcall = './.' if ($sampleAlleleCalls->[0] =~ /\S\S/);
148 }
149 $info_NS-- if ($sampleGTcall =~ /\./);
150 push (@allGTcalls, $sampleGTcall);
151 $info_AF[$1]++ if ($sampleGTcall =~ /^(\d+)/);
152 $info_AF[$1]++ if ($sampleGTcall =~ /^\d+\/(\d+)$/);
153 }
154 $filter = 's50' if ($info_NS/scalar @$samples < 0.5);
155 print "$filter\t";
156 my ($numAlleles) = 0;
157 foreach (@info_AF) {$numAlleles += $_;}
158 $numAlleles = 1 if ($numAlleles == 0); #In some loci, all samples are no-called
159 my (@alleleFreq) = map (sprintf ('%.2f', $_/$numAlleles), @info_AF);
160 shift @alleleFreq; #don't need to print the reference allele frequency
161 print "NS=$info_NS;RSID=", join (",", @info_RSid), ";AF=", join (",", @alleleFreq), ";$info_DB\t$format\t", join ("\t", @allGTcalls), "\n";
162 }
163
164 #Convert a list of allele calls to diploid VCF genotype calls
165 #[11] becomes 1/1
166 #[01] becomes 0/1
167 #[00] becomes 0/0
168 #[11, 00] becomes 1/1
169 #[01, 01] becomes 1/2
170 #[00, 01, 1N] becomes 2/3
171 #[01, 01, 01] becomes ./. (unknown)
172 #[1N, 00, 00, 00] becomes 1/.
173 sub mkDiploidGTcall {
174 my ($sampleAlleleCalls) = @_;
175 #[11, 00, NN, 00], [1N, NN, 01, 00], [1N, NN, 1N, 00], [01, NN, 01, 00], [1N, 00, 00, 00]
176 my ($numOfOnes) = 0;
177 my ($sampleGTcall) = './.'; #diploid no-call in VCF
178 my (%merged);
179 for (my $i=0; $i<scalar @$sampleAlleleCalls; $i++) {
180 push (@{$merged{$sampleAlleleCalls->[$i]}}, $i+1);
181 my ($ones) = $sampleAlleleCalls->[$i] =~ tr/1/1/;
182 $numOfOnes += $ones;
183 }
184 return 'warning-more_than_two_1s' if ($numOfOnes > 2);
185 return $merged{'11'}->[0] . '/' . $merged{'11'}->[0] if (exists $merged{'11'});
186 if (exists $merged{'01'}) {
187 return $merged{'01'}->[0] . '/' . $merged{'01'}->[1] if (defined $merged{'01'}->[1]);
188 return $merged{'01'}->[0] . '/' . $merged{'1N'}->[0] if (exists $merged{'1N'});
189 return '0/' . $merged{'01'}->[0] if (!exists $merged{'NN'} && !exists $merged{'0N'});
190 return $merged{'01'}->[0] . '/.';
191 }
192 if (exists $merged{'1N'}) {
193 return $merged{'1N'}->[0] . '/' . $merged{'1N'}->[1] if (defined $merged{'1N'}->[1]);
194 return $merged{'1N'}->[0] . '/.';
195 }
196 return '0/0' if (!exists $merged{'NN'} && !exists $merged{'0N'});
197 return '0/.' if (exists $merged{'0N'} && scalar @{$merged{'0N'}} == 1);
198 return $sampleGTcall;
199 }
200
201 #Convert a list of allele calls to haploid VCF genotype calls
202 #[1] becomes 1
203 #[N] becomes .
204 #[0] becomes 0
205 #[1, 0] becomes 1
206 #[0, 1, N] becomes 2
207 #[1, 1] becomes .
208 sub mkHaploidGTcall {
209 my ($sampleAlleleCalls) = @_;
210 my ($sampleGTcall) = '.'; #haploid no-call in VCF
211 my (%merged);
212 for (my $i=0; $i<scalar @$sampleAlleleCalls; $i++) {
213 push (@{$merged{$sampleAlleleCalls->[$i]}}, $i+1);
214 }
215 if (exists $merged{'1'}) {
216 return "warning-more_than_one_1_haploid" if (defined $merged{'1'}->[1]);
217 return $merged{'1'}->[0];
218 }
219 return '0' if (!exists $merged{'N'});
220 return $sampleGTcall;
221 }