Mercurial > repos > miller-lab > snp_analysis_conversion

comparison snp_analysis_conversion/master2gd_snp.pl @ 0:3871157bc013

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial upload to toolshed.g2 via UI.
author cathy
date Tue, 28 May 2013 17:01:14 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:3871157bc013
1 #!/usr/bin/perl -w
2 use strict;
3
4 #convert from master variant file to snp table (Webb format?)
5 #new format for version 2.0, also different format for cancer normal pairs
6 #set columns for 2.0 version Cancer format
7 my $aCnt1 = 21;
8 my $aCnt2 = 22;
9
10 #snp table format:
11 #1. chr
12 #2. position (0 based)
13 #3. ref allele
14 #4. second allele
15 #5. overall quality
16 #foreach individual (6-9, 10-13, ...)
17 #a. count of allele in 3
18 #b. count of allele in 4
19 #c. genotype call (-1, or count of ref allele)
20 #d. quality of genotype call (quality of non-ref allele from masterVar)
21
22 if (!@ARGV) {
23 print "usage: master2gd_snp.pl masterVar.txt[.gz|.bz2] [-tab=snpTable.txt -addColsOnly -build=hg19 -name=na ] > newSnpTable.txt\n";
24 exit;
25 }
26
27 my $in = shift @ARGV;
28 my $tab;
29 my $tabOnly;
30 my $build;
31 my $name;
32 foreach (@ARGV) {
33 if (/-tab=(.*)/) { $tab = $1; }
34 elsif (/-addColsOnly/) { $tabOnly = 1; }
35 elsif (/-build=(.*)/) { $build = $1; }
36 elsif (/-name=(.*)/) { $name = $1; }
37 }
38
39 #WARNING loads snp table in memory, this could take > 1G ram
40 my %old;
41 my $colcnt = 0;
42 my @head;
43 if ($tab) {
44 open(FH, $tab) or die "Couldn't open $tab, $!\n";
45 while (<FH>) {
46 chomp;
47 if (/^#/) { push(@head, $_); next; }
48 my @f = split(/\t/);
49 $old{"$f[0]:$f[1]"} = join("\t", @f);
50 $colcnt = scalar @f;
51 }
52 close FH or die "Couldn't close $tab, $!\n";
53 }
54
55 if ($in =~ /.gz$/) {
56 open(FH, "zcat $in |") or die "Couldn't open $in, $!\n";
57 }elsif ($in =~ /.bz2$/) {
58 open(FH, "bzcat $in |") or die "Couldn't open $in, $!\n";
59 }else {
60 open(FH, $in) or die "Couldn't open $in, $!\n";
61 }
62 prepHeader();
63 if (@head) { #keep old header, add new?
64 print join("\n", @head), "\n";
65 }
66 while (<FH>) {
67 chomp;
68 #FORMAT_VERSION 2.0
69 if (/^#FORMAT_VERSION\s+1\./) {
70 $aCnt1 = 16;
71 $aCnt2 = 17;
72 }
73 if (/^#/) { next; }
74 if (/^>/) { next; } #headers
75 if (/^\s*$/) { next; }
76 my @f = split(/\t/);
77 if (!$f[6]) { next; } #WHAT? most likely still zipped?
78 if ($f[6] ne 'snp') { next; } #table only has substitutions
79 if ($f[5] eq 'het-alt') { next; } #skip heterozygous with no ref match
80 if ($f[5] =~ /(hom|het)/) { #zygosity #haploid chrX and chrY?
81 my $a = $f[7]; #reference allele
82 my $a2;
83 my $freq;
84 my $freq2;
85 my $sc = -1;
86 my $alt;
87 my $g = 1; #genotype == ref allele count
88 if ($f[8] eq $f[9]) { #should be homozygous?
89 $a2 = $f[8];
90 $g = 0;
91 if ($f[10] && $f[10] ne '') { $sc = $f[10]; }#is this the best one to use? or smallest?
92 }else {
93 if ($a ne $f[8]) { $a2 = $f[8]; $alt = 8; }
94 elsif ($a ne $f[9]) { $a2 = $f[9]; $alt = 9; }
95 }
96 if (defined $f[10] && defined $f[11] && $alt) { #VAF score in 2.0 format
97 if ($f[$alt+2] && $f[$alt+2] ne '') { $sc = $f[$alt+2]; }
98 }
99 #version 1.12 columns 16 & 17, version 2.0 Cancer columns 21 & 22
100 if (defined $f[$aCnt1] && defined $f[$aCnt2] && $alt) {
101 if ($alt == 8) {
102 $freq = $f[$aCnt2];
103 $freq2 = $f[$aCnt1];
104 }elsif ($alt == 9) {
105 $freq = $f[$aCnt1];
106 $freq2 = $f[$aCnt2];
107 }
108 }elsif (defined $f[$aCnt1]) {
109 $freq = 0;
110 $freq2 = $f[$aCnt1];
111 }
112 #if starting a new table or new SNP in old table
113 #add option to only build on current table?
114 if (!$tab) {
115 print "$f[2]\t$f[3]\t$a\t$a2\t-1";
116 }elsif (!$tabOnly && !exists $old{"$f[2]:$f[3]"}) {
117 print "$f[2]\t$f[3]\t$a\t$a2\t-1";
118 }elsif (exists $old{"$f[2]:$f[3]"}) {
119 print $old{"$f[2]:$f[3]"};
120 $old{"$f[2]:$f[3]"} = ''; #unset so we know it is printed
121 }elsif ($tabOnly && !exists $old{"$f[2]:$f[3]"}) {
122 next; #skip this one entirely
123 }
124 if ($colcnt && !exists $old{"$f[2]:$f[3]"}) {
125 #new SNP pad for missing individuals
126 my $i = 5;
127 while ($i < $colcnt) {
128 print "\t-1\t-1\t-1\t-1";
129 $i += 4;
130 }
131 }
132 #add columns for individual
133 print "\t$freq\t$freq2\t$g\t$sc\n";
134 }elsif ($f[5] eq 'hap') {
135 my $g = 0;
136 my $freq = 0;
137 my $freq2 = 0;
138 if (defined $f[10]) { $freq2 = $f[10]; }
139 my $sc = -1;
140 if (defined $f[$aCnt1] && $f[$aCnt1] ne '') { $sc = $f[$aCnt1]; }
141 if ($f[8]) {
142 if (!$tab) {
143 print "$f[2]\t$f[3]\t$f[7]\t$f[8]\t-1";
144 }elsif (!$tabOnly && !exists $old{"$f[2]:$f[3]"}) {
145 print "$f[2]\t$f[3]\t$f[7]\t$f[8]\t-1";
146 }elsif (exists $old{"$f[2]:$f[3]"}) {
147 print $old{"$f[2]:$f[3]"};
148 $old{"$f[2]:$f[3]"} = ''; #unset so we know it is printed
149 }elsif ($tabOnly && !exists $old{"$f[2]:$f[3]"}) {
150 next; #skip this one entirely
151 }
152 if ($colcnt && !exists $old{"$f[2]:$f[3]"}) {
153 #new SNP pad for missing individuals
154 my $i = 5;
155 while ($i < $colcnt) {
156 print "\t-1\t-1\t-1\t-1";
157 $i += 4;
158 }
159 }
160 #add columns for individual
161 print "\t$freq\t$freq2\t$g\t$sc\n";
162 }
163 }
164 }
165 close FH or die "Couldn't close $in, $!\n";
166
167 #if adding to a snp table, now we need to finish those not in the latest set
168 foreach my $k (keys %old) {
169 if ($old{$k} ne '') { #not printed yet
170 print $old{$k}, "\t-1\t-1\t-1\t-1\n"; #plus blank for this one
171 }
172 }
173
174 exit;
175
176 #parse old header and add or create new
177 sub prepHeader {
178 if (!$build) { $build = 'hg19'; } #set default
179 my @cnames;
180 my @ind;
181 my $n;
182 if (@head) { #parse previous header
183 my $h = join("", @head); #may split between lines
184 if ($h =~ /"column_names":\[(.*?)\]/) {
185 my @t = split(/,/, $1);
186 foreach (@t) { s/"//g; }
187 @cnames = @t;
188 $n = $cnames[$#cnames];
189 $n =~ s/Q//;
190 $n++;
191 }
192 if ($h =~ /"dbkey":"(.*?)"/) { $build = $1; }
193 if ($h =~ /"individuals":\[(.*)\]/) {
194 my $t = $1;
195 $t =~ s/\]\].*/]/; #remove if there is more categories
196 @ind = split(/,/, $t);
197 }
198 }else { #start new header
199 @cnames = ("chr", "pos", "A", "B", "Q");
200 $n = 1;
201 }
202 #add current
203 if (!$name) { $name= 'na'; }
204 my $stcol = $colcnt + 1;
205 if ($stcol == 1) { $stcol = 6; } #move past initial columns
206 push(@ind, "[\"$name\",$stcol]");
207 push(@cnames, "${n}A", "${n}B", "${n}G", "${n}Q");
208 #reassign head
209 undef @head;
210 foreach (@cnames) { $_ = "\"$_\""; } #quote name
211 $head[0] = "#{\"column_names\":[" . join(",", @cnames) . "],";
212 $head[1] = "#\"individuals\":[" . join(",", @ind) . "],";
213 $head[2] = "#\"dbkey\":\"$build\",\"pos\":2,\"rPos\":2,\"ref\":1,\"scaffold\":1,\"species\":\"$build\"}";
214 }
215 ####End
216
217 ##example header
218 #{"column_names":["chr","pos","A","B","Q","1A","1B","1G","1Q","2A","2B","2G","2Q","3A","3B","3G","3Q","4A","4B","4G","4Q","5A","5B","5G","5Q","6A","6B","6G","6Q","7A","7B","7G","7Q","8A","8B","8G",
219 #"8Q","9A","9B","9G","9Q","10A","10B","10G","10Q"],"dbkey":"hg19","individuals":[["Boh_15M",6],["Boh_19M",10],["Paya_27F",14],["Paya_2F",18],["Paya_32F",22],["Ruil_2M",26],["Ruil_36M",30],["Ruil_3M",
220 #34],["Ruil_40",38],["Ruil_47F",42]],"pos":2,"rPos":2,"ref":1,"scaffold":1,"species":"hg19"}
221 #chr1 10290 C T 46.4 0 2 0 7 1 2 0 4 3 2 1 22 0 0 -1 0 1 0 1 4 0 2 0 7 0 0 -1 0 2 3 1 14 0 1 0 4 1 1 1 6