0
|
1 #!/usr/bin/perl -w
|
|
2 use strict;
|
|
3
|
|
4 #convert from gd_snp file to vcf file (with dbSNP fields)
|
|
5
|
|
6 #gd_snp table format:
|
|
7 #1. chr
|
|
8 #2. position (0 based)
|
|
9 #3. ref allele
|
|
10 #4. second allele
|
|
11 #5. overall quality
|
|
12 #foreach individual (6-9, 10-13, ...)
|
|
13 #a. count of allele in 3
|
|
14 #b. count of allele in 4
|
|
15 #c. genotype call (-1, or count of ref allele)
|
|
16 #d. quality of genotype call (quality of non-ref allele from masterVar)
|
|
17
|
|
18 if (!@ARGV) {
|
|
19 print "usage: gd_snp2vcf.pl file.gd_snp[.gz|.bz2] -geno=8[,12:16,20...] -handle=HANDLE -batch=BATCHNAME -ref=REFERENCEID [-bioproj=XYZ -biosamp=ABC -pop=POPID[,POPID2...] -chrCol=9 -posCol=9 ] > snpsForSubmission.vcf\n";
|
|
20 exit;
|
|
21 }
|
|
22
|
|
23 my $in = shift @ARGV;
|
|
24 my $genoCols = '';
|
|
25 my $handle;
|
|
26 my $batch;
|
|
27 my $bioproj;
|
|
28 my $biosamp;
|
|
29 my $ref;
|
|
30 my $pop;
|
|
31 my $cr = 0; #allow to use alternate reference?
|
|
32 my $cp = 1;
|
|
33 my $meta;
|
|
34 foreach (@ARGV) {
|
|
35 if (/-geno=([0-9,]+)/) { $genoCols .= "$1:"; }
|
|
36 elsif (/-handle=(.*)/) { $handle = $1; }
|
|
37 elsif (/-batch=(.*)/) { $batch = $1; }
|
|
38 elsif (/-bioproj=(.*)/) { $bioproj = $1; }
|
|
39 elsif (/-biosamp=(.*)/) { $biosamp = $1; }
|
|
40 elsif (/-ref=(.*)/) { $ref = $1; }
|
|
41 elsif (/-population=(\S+)/) { $pop = $1; }
|
|
42 elsif (/-chrCol=(\d+)/) { $cr = $1 - 1; }
|
|
43 elsif (/-posCol=(\d+)/) { $cp = $1 - 1; }
|
|
44 elsif (/-metaOut=(.*)/) { $meta = $1; }
|
|
45 }
|
|
46 if ($cr < 0 or $cp < 0) { die "ERROR the column numbers should be 1 based.\n"; }
|
|
47
|
|
48 #remove trailing delimiters
|
|
49 $genoCols =~ s/,:/:/g;
|
|
50 $genoCols =~ s/[,:]$//;
|
|
51
|
|
52 my @gnc = split(/,|:/, $genoCols);
|
|
53
|
|
54 if ($in =~ /.gz$/) {
|
|
55 open(FH, "zcat $in |") or die "Couldn't open $in, $!\n";
|
|
56 }elsif ($in =~ /.bz2$/) {
|
|
57 open(FH, "bzcat $in |") or die "Couldn't open $in, $!\n";
|
|
58 }else {
|
|
59 open(FH, $in) or die "Couldn't open $in, $!\n";
|
|
60 }
|
|
61 my @head = prepHeader();
|
|
62 if (@head) {
|
|
63 print join("\n", @head), "\n";
|
|
64 #now column headers
|
|
65 print "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO";
|
|
66 if (defined $pop) {
|
|
67 $pop =~ s/,$//;
|
|
68 my $t = $pop;
|
|
69 $t =~ s/,/\t/g;
|
|
70 print "\tFORMAT\t$t";
|
|
71 }
|
|
72 print "\n";
|
|
73 }
|
|
74 while (<FH>) {
|
|
75 chomp;
|
|
76 if (/^#/) { next; }
|
|
77 if (/^\s*$/) { next; }
|
|
78 my @f = split(/\t/);
|
|
79 #vcf columns: chrom pos id ref alt qual filter info
|
|
80 # info must have VRT=[0-9] 1==SNV 2=indel 6=NoVariation 8=MNV ...
|
|
81 my $vrt = 1;
|
|
82 if ($f[2] !~ /^[ACTG]$/ or $f[3] !~ /^[ACTG]$/) {
|
|
83 die "Sorry this can only do SNV's at this time\n";
|
|
84 }
|
|
85 if (scalar @gnc == 1) { #single genotype column
|
|
86 if (!defined $f[4] or $f[4] == -1) { $f[4] = '.'; }
|
|
87 if ($f[$gnc[0]-1] == 2) { $vrt = 6; } #reference match
|
|
88 print "$f[$cr]\t$f[$cp]\t$f[$cr];$f[$cp]\t$f[2]\t$f[3]\t$f[4]\t.\tVRT=$vrt\n";
|
|
89 #TODO? put read counts in comment?
|
|
90 }elsif ($pop) { #do as population
|
|
91 my @cols;
|
|
92 foreach my $gp (split(/:/,$genoCols)) { #foreach population
|
|
93 my @g = split(/,/, $gp);
|
|
94 my $totChrom = 2*(scalar @g);
|
|
95 my $totRef = 0;
|
|
96 foreach my $i (@g) { if ($f[$i-1] == -1) { next; } $totRef += $f[$i-1]; }
|
|
97 if ($totChrom == $totRef) { $vrt = 6; }
|
|
98 if ($totRef > $totChrom) { die "ERROR likely the wrong column was chosen for genotype\n"; }
|
|
99 my $altCnt = $totChrom - $totRef;
|
|
100 push(@cols, "$totChrom:$altCnt");
|
|
101 }
|
|
102 print "$f[$cr]\t$f[$cp]\t$f[$cr];$f[$cp]\t$f[2]\t$f[3]\t$f[4]\t.\tVRT=$vrt\tNA:AC\t", join("\t", @cols), "\n";
|
|
103 }else { #leave allele counts off
|
|
104 my $totChrom = 2*(scalar @gnc);
|
|
105 my $totRef = 0;
|
|
106 foreach my $i (@gnc) { if ($f[$i-1] == -1) { next; } $totRef += $f[$i-1]; }
|
|
107 if ($totChrom == $totRef) { $vrt = 6; }
|
|
108 print "$f[$cr]\t$f[$cp]\t$f[$cr];$f[$cp]\t$f[2]\t$f[3]\t$f[4]\t.\tVRT=$vrt\n";
|
|
109 }
|
|
110 }
|
|
111 close FH or die "Couldn't close $in, $!\n";
|
|
112
|
|
113 if ($meta) {
|
|
114 open(FH, ">", $meta) or die "Couldn't open $meta, $!\n";
|
|
115 print FH "TYPE: CONT\n",
|
|
116 "HANDLE: $handle\n",
|
|
117 "NAME: \n",
|
|
118 "FAX: \n",
|
|
119 "TEL: \n",
|
|
120 "EMAIL: \n",
|
|
121 "LAB: \n",
|
|
122 "INST: \n",
|
|
123 "ADDR: \n",
|
|
124 "||\n",
|
|
125 "TYPE: METHOD\n",
|
|
126 "HANDLE: $handle\n",
|
|
127 "ID: \n",
|
|
128 "METHOD_CLASS: Sequence\n",
|
|
129 "TEMPLATE_TYPE: \n",
|
|
130 "METHOD:\n",
|
|
131 "||\n";
|
|
132 if ($pop) {
|
|
133 my @p = split(/,/, $pop);
|
|
134 foreach my $t (@p) {
|
|
135 print FH
|
|
136 "TYPE: POPULATION\n",
|
|
137 "HANDLE: $handle\n",
|
|
138 "ID: $t\n",
|
|
139 "POPULATION: \n",
|
|
140 "||\n";
|
|
141 }
|
|
142 }
|
|
143 print FH "TYPE: SNPASSAY\n",
|
|
144 "HANDLE: $handle\n",
|
|
145 "BATCH: $batch\n",
|
|
146 "MOLTYPE: \n",
|
|
147 "METHOD: \n",
|
|
148 "ORGANISM: \n",
|
|
149 "||\n",
|
|
150 "TYPE: SNPPOPUSE | SNPINDUSE\n",
|
|
151 "HANDLE: $handle\n",
|
|
152 "BATCH: \n",
|
|
153 "METHOD: \n",
|
|
154 "||\n";
|
|
155
|
|
156 close FH or die "Couldn't close $meta, $!\n";
|
|
157 }
|
|
158
|
|
159 exit 0;
|
|
160
|
|
161 #parse old header and add or create new
|
|
162 sub prepHeader {
|
|
163 my @h;
|
|
164 $h[0] = '##fileformat=VCFv4.1';
|
|
165 my ($day, $mo, $yr) = (localtime)[3,4,5];
|
|
166 $mo++;
|
|
167 $yr+=1900;
|
|
168 $h[1] = '##fileDate=' . "$yr$mo$day";
|
|
169 $h[2] = "##handle=$handle";
|
|
170 $h[3] = "##batch=$batch";
|
|
171 my $i = 4;
|
|
172 if ($bioproj) { $h[$i] = "##bioproject_id=$bioproj"; $i++; }
|
|
173 if ($biosamp) { $h[$i] = "##biosample_id=$biosamp"; $i++; }
|
|
174 $h[$i] = "##reference=$ref"; ##reference=GCF_999999.99
|
|
175 #$i++;
|
|
176 #$h[$i] = '##INFO=<ID=LID, Number=1,Type=string, Description="Unique local variation ID or name for display. The LID provided here combined with the handle must be unique for a particular submitter.">'
|
|
177 $i++;
|
|
178 $h[$i] = '##INFO=<ID=VRT,Number=1,Type=Integer,Description="Variation type,1 - SNV: single nucleotide variation,2 - DIV: deletion/insertion variation,3 - HETEROZYGOUS: variable, but undefined at nucleotide level,4 - STR: short tandem repeat (microsatellite) variation, 5 - NAMED: insertion/deletion variation of named repetitive element,6 - NO VARIATON: sequence scanned for variation, but none observed,7 - MIXED: cluster contains submissions from 2 or more allelic classes (not used) ,8 - MNV: multiple nucleotide variation with all eles of common length greater than 1,9 - Exception">';
|
|
179 #sometimes have allele freqs?
|
|
180 if (defined $pop) {
|
|
181 $i++;
|
|
182 $h[$i] = "##FORMAT=<ID=NA,Number=1,Type=Integer,Description=\"Number of alleles for the population.\"";
|
|
183 $i++;
|
|
184 $h[$i] = '##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count for each alternate allele.">';
|
|
185 my @p = split(/,/, $pop);
|
|
186 foreach my $t (@p) {
|
|
187 $i++;
|
|
188 $h[$i] = "##population_id=$t";
|
|
189 }
|
|
190 }
|
|
191 #PMID?
|
|
192 ##INFO=<ID=PMID,Number=.,Type=Integer,Description="PubMed ID linked to variation if available.">
|
|
193
|
|
194 return @h;
|
|
195 }
|
|
196 ####End
|
|
197
|