bin/convert_vcf_to_hmp.pl

   1
   2 use strict;
   3 use Data::Dumper;
   4
   5 my $vcf_file = shift;
   6
   7 #print STDERR Dumper($vcf_file);
   8
   9 open(my $F, "<", $vcf_file) || die "Can't find file $vcf_file";
  10
  11 open(my $G, ">", $vcf_file.".hmp") || die "Can't open $vcf_file.hmp";
  12
  13 my $line;
  14 while (<$F>) {
  15     my $line = $_;
  16     if ($line =~ m/^\#\#/) {
  17         print STDERR "SKIPPING ## lines...\n";
  18     }
  19     else {
  20         last;
  21     }
  22 }
  23
  24 my $header = $line;
  25
  26 print STDERR $header;
  27 chomp($header);
  28
  29 my @keys = split("\t", $header);
  30 #print STDERR Dumper($keys[1]);
  31
  32 for(my $n=0; $n <@keys; $n++) {
  33     if ($keys[$n] =~ /\|CO\_/) {
  34         $keys[$n] =~ s/\|CO\_.*//;
  35     }
  36 }
  37 my @data = ();
  38
  39 my %nuconv = (
  40     'A/A' => 'A',
  41     'G/G' => 'G',
  42     'T/T' => 'T',
  43     'C/C' => 'C',
  44     'A/G' => 'R',
  45     'G/A' => 'R',
  46     'C/T' => 'Y',
  47     'T/C' => 'Y',
  48     'G/C' => 'S',
  49     'C/G' => 'S',
  50     'A/T' => 'W',
  51     'T/A' => 'W',
  52     'G/T' => 'K',
  53     'T/G' => 'K',
  54     'A/C' => 'M',
  55     'C/A' => 'M',
  56     './.' => 'N',
  57     );
  58
  59
  60 while (<$F>) {
  61     chomp;
  62     my %line;
  63     my @fields = split /\t/;
  64
  65     for(my $n=0; $n <@keys; $n++) {
  66         if (exists($fields[$n]) && defined($fields[$n])) {
  67             $line{$keys[$n]}=$fields[$n];
  68         }
  69     }
  70     push @data, \%line;
  71 }
  72
  73 foreach my $line (@data) {
  74     my @formats = split /\:/, $line->{FORMAT};
  75     my $gtindex = 0;
  76     for(my $n =0; $n<@formats; $n++) {
  77         if ($formats[$n] eq 'GT') {
  78             $gtindex=$n;
  79         }
  80     }
  81
  82     for(my $gt = 0; $gt < @keys; $gt++) {
  83         my @scores = split /\:/, $line->{$gt};
  84         my $genotype;
  85         if ($scores[$gtindex] eq '0/0') {
  86             $genotype = $line->{REF}."/".$line->{REF};
  87         }
  88         if ($scores[$gtindex] eq '0/1') {
  89             $genotype = $line->{REF}."/".$line->{ALT};
  90         }
  91         if ($scores[$gtindex] eq '1/1') {
  92             $genotype = $line->{ALT}."/".$line->{ALT};
  93         }
  94
  95         my $genotype_hmp = $nuconv{$genotype};
  96     }
  97 }
  98
  99
 100 close($F);
 101 close($G);