lib/CXGN/Genotype/SNP.pm

   1
   2 package CXGN::Genotype::SNP;
   3
   4 use Moose;
   5
   6 use Math::BigInt;
   7
   8 sub BUILD {
   9     my $self = shift;
  10     my $args = shift;
  11
  12     if ($args->{vcf_string}) {
  13         #print STDERR "Building SNP from vcf_string $args->{vcf_string}...\n";
  14         $self->from_vcf_string($args->{vcf_string});
  15
  16         #print STDERR "Counts: ".$self->ref_count().", ".$self->alt_count()."\n";
  17     }
  18 }
  19
  20 has 'id' => (isa => 'Str',
  21              is => 'rw',
  22     );
  23
  24 has 'format' => ( isa => 'Str',
  25                   is => 'rw',
  26     );
  27 has 'vcf_string' => ( isa => 'Str',
  28                       is  => 'rw',
  29     );
  30 has 'accession' => (isa => 'Str',
  31                     is => 'rw',
  32     );
  33
  34 has 'ref_allele' => ( isa => 'Str',
  35                    is  => 'rw',
  36     );
  37
  38 has 'alt_allele' => ( isa => 'Str',
  39                    is  => 'rw',
  40     );
  41
  42 has 'ref_count' => ( isa => 'Int',
  43                   is  => 'rw',
  44     );
  45
  46 has 'alt_count' => ( isa => 'Int',
  47                   is  => 'rw',
  48     );
  49
  50 has 'dosage' => (
  51                   is  => 'rw',
  52     );
  53
  54
  55 sub from_vcf_string {
  56     my $self = shift;
  57     my $raw = shift;
  58
  59     my @ids;
  60     if ($self->format()) {
  61         @ids = split ":", $self->format();
  62     }
  63     if (!@ids) {
  64         # usual order... but this is just guesswork
  65         #print STDERR "(Guessing the format order...)\n";
  66         @ids = ( 'GT', 'AD', 'DP', 'GQ', 'PL' );
  67     }
  68
  69     my @values = split /\:/, $raw;
  70
  71     my %fields;
  72     for(my $i=0; $i<@ids; $i++) {
  73          $fields{$ids[$i]} = $values[$i];
  74      }
  75
  76      #my ($allele, $counts) = split /\:/, $raw;
  77      my ($a1, $a2) = ("","");
  78      if (!exists($fields{GT})) {
  79          print STDERR "No allele calls found for snp ".$self->id()."\n";
  80      }
  81      else {
  82          ($a1, $a2) = split /\//, $fields{GT};
  83      }
  84      $self->ref_allele($a1);
  85      $self->alt_allele($a2);
  86
  87      my ($c1, $c2);
  88      if (!exists($fields{AD})) {
  89          $c1 = 0;
  90          #print STDERR "C1: $c1\n";
  91          $c2 = 0;
  92          #print STDERR "C2: $c2\n";
  93      }
  94      else {
  95          my @counts = split /\,/, $fields{AD};
  96          # print STDERR "FIELDS: $fields{AD}\n";
  97          if (@counts > 2) {
  98              #print STDERR "Multiple alleles found for SNP ".$self->id()."\n";
  99          }
 100          ($c1, $c2) = @counts;
 101      }
 102     if (!defined($c1)) { $c1=0; }
 103     if (!defined($c2)) { $c2=0; }
 104      $self->ref_count($c1);
 105      $self->alt_count($c2);
 106
 107     # debug
 108     #if ($self->id() eq "1002:250060174") {
 109         #print STDERR $self->id().": ".$fields{GT}.", ".$fields{AD}."\n";
 110     #}
 111
 112     return ($c1, $c2);
 113 }
 114
 115 sub good_call {
 116     my $self = shift;
 117     my $call_sum_min = shift || 2;
 118     my ($c1, $c2) = ($self->ref_count(), $self->alt_count());
 119     if ( ($c1 + $c2) < $call_sum_min) {
 120         return 0;
 121     }
 122     return 1;
 123 }
 124
 125
 126 __PACKAGE__->meta->make_immutable;
 127
 128 1;