Mercurial > repos > miller-lab > snp_analysis_conversion

comparison vcf2pgSnpMult.pl @ 2:35c20b109be5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Retrying upload with "bare" tarball (i.e. one without a top containing directory).
author cathy
date Tue, 28 May 2013 17:54:02 -0400
parents
children
comparison
equal deleted inserted replaced
1:1d8b23a21735 2:35c20b109be5
1 #!/usr/bin/perl -w
2 use strict;
3
4 #convert from a vcf file to a pgSnp file with multiple sets of the allele
5 # specific columns
6 #frequency count = chromosome count
7
8 my $in;
9 my $stCol = 9;
10 my $endCol;
11 if (@ARGV && scalar @ARGV == 1) {
12 $in = shift @ARGV;
13 }else {
14 print "usage: vcf2pgSnpMult.pl file.vcf > file.pgSnpMult\n";
15 exit;
16 }
17
18 if ($in =~ /.gz$/) {
19 open(FH, "zcat $in |") or die "Couldn't open $in, $!\n";
20 }else {
21 open(FH, $in) or die "Couldn't open $in, $!\n";
22 }
23 while (<FH>) {
24 chomp;
25 if (/^\s*#/) { next; } #skip comments/headers
26 if (/^\s*$/) { next; } #skip blank lines
27 my @f = split(/\t/);
28 #chr pos1base ID refNt altNt[,|D#|Int] quality filter info format geno1 ...
29 my $a;
30 my %nt;
31 my %all;
32 my $cnt = 0;
33 my $var;
34 if ($f[3] eq 'N') { next; } #ignore ref=N
35 if ($f[4] =~ /[DI]/ or $f[3] =~ /[DI]/) { next; } #don't do microsatellite
36 if ($f[6] && !($f[6] eq '.' or $f[6] eq 'PASS')) { next; } #filtered for some reason
37 my $ind = 0;
38 if ($f[8] ne 'GT') { #more than just genotype
39 my @t = split(/:/, $f[8]);
40 foreach (@t) { if ($_ eq 'GT') { last; } $ind++; }
41 if ($ind == 0 && $f[8] !~ /^GT/) { die "ERROR couldn't find genotype in format $f[8]\n"; }
42 }
43 if (!$endCol) { $endCol = $#f; }
44 #put f[3] => nt{0} and split f[4] for rest of nt{}
45 $nt{0} = $f[3];
46 my @t = split(/,/, $f[4]);
47 for (my $i=0; $i<=$#t; $i++) {
48 my $j = $i + 1;
49 $nt{$j} = $t[$i];
50 }
51 if ($f[0] !~ /chr/) { $f[0] = "chr$f[0]"; }
52 print "$f[0]\t", ($f[1]-1), "\t$f[1]"; #position info
53 foreach my $col ($stCol .. $endCol) { #add each individual (4 columns)
54 if ($ind > 0) {
55 my @t = split(/:/, $f[$col]);
56 $f[$col] = $t[$ind] . ":"; #only keep genotype part
57 }
58 print "\t";
59 if ($f[$col] =~ /^(\d).(\d)/) {
60 my $a1 = $1;
61 my $a2 = $2;
62 if (!exists $nt{$a1}) { die "ERROR bad allele $a1 in $f[3] $f[4]\n"; }
63 if (!exists $nt{$a2}) { die "ERROR bad allele $a2 in $f[3] $f[4]\n"; }
64 if ($a1 eq $a2) { #homozygous
65 print "$nt{$a1}\t1\t2\t0";
66 }else { #heterozygous
67 print "$nt{$a1}/$nt{$a2}\t2\t1,1\t0,0";
68 }
69 }elsif ($f[$col] =~ /^(\d):/) { #chrY or male chrX, single
70 my $a1 = $1;
71 if (!exists $nt{$a1}) { die "ERROR bad allele $a1 in $f[3] $f[4]\n"; }
72 print "$nt{$a1}\t1\t1\t0";
73 }else { #don't know how to parse
74 die "ERROR unknown genotype $f[$col]\n";
75 }
76 }
77 print "\n"; #end this SNP
78 }
79 close FH;
80
81 exit;