lib/CXGN/SNPs.pm

   1
   2 =head1 NAME
   3
   4 CXGN::SNPs - deal with collections of SNPs
   5
   6 =head1 SYNOPSYS
   7
   8  my $snps = CXGN::SNPs->new( { file => $file });
   9  my $allele_feq = $snps->allele_freq();
  10  my $dosage = $snps->calculate_snp_dosage($snp); # provide a CXGN::Genotype::SNP object
  11  # etc
  12
  13 =head1 AUTHOR
  14
  15 Lukas Mueller <lam87@cornell.edu>
  16
  17 =head1 METHODS
  18
  19 =cut
  20
  21 package CXGN::SNPs;
  22
  23 use Moose;
  24
  25 use Data::Dumper;
  26 use Math::BigInt;
  27
  28 =head2 accessors id()
  29
  30 The id of the SNP.
  31
  32 =cut
  33
  34 has 'id' => ( isa => 'Str',
  35               is => 'rw',
  36     );
  37
  38 =head2 accessors accessions()
  39
  40 All the accessions parsed in the SNP
  41
  42 =cut
  43
  44 has 'accessions' => ( isa => 'ArrayRef',
  45                        is  => 'rw',
  46     );
  47
  48 =head2 accessors depth()
  49
  50 The total depth of the SNP (usually from VCF file)
  51
  52 =cut
  53
  54 has 'depth' => ( isa => 'Int',
  55                  is => 'rw',
  56     );
  57
  58 =head2 accessors ignore_accessions()
  59
  60 A hashref with accessions that have not passed some quality check and that should be ignored.
  61
  62 =cut
  63
  64 has 'ignore_accessions' => (isa => 'HashRef',
  65                             is => 'rw',
  66     );
  67
  68 =head2 accessors valid_accessions()
  69
  70 A ArrayRef of valid accessions (the complement set to ignore_accessions)
  71
  72 =cut
  73
  74 has 'valid_accessions' => (isa => 'ArrayRef',
  75                            is => 'rw',
  76     );
  77
  78
  79 has 'scores'  => ( isa => 'HashRef',
  80                        is  => 'rw',
  81     );
  82
  83 #has 'bad_clones'  => ( isa => 'ArrayRef',
  84 #                      is  => 'rw',
  85  #   );
  86
  87 =head2 accessors snps()
  88
  89 An HashRef of CXGN::Genotype::SNP objects that represent the SNP calls, with genotype ids as the hash key.
  90
  91 =cut
  92
  93 has 'snps' => ( isa => 'HashRef',
  94                 is  => 'rw',
  95     );
  96
  97 =head2 accessors maf()
  98
  99 The minor allele frequency of the SNP set.
 100
 101 =cut
 102
 103 has 'maf'  => ( isa => 'Num',
 104                 is  => 'rw',
 105                 default => sub { 0.999 },
 106     );
 107
 108 =head2 accessors allele_freq()
 109
 110 The allele frequency of the SNP set.
 111
 112 =cut
 113
 114 has 'allele_freq' => ( isa => 'Num',
 115                        is  => 'rw',
 116     );
 117
 118
 119
 120 has 'chr' => ( isa => 'Str',
 121                is => 'rw',
 122     );
 123
 124 has 'position' => (isa => 'Int',
 125                    is => 'rw',
 126     );
 127
 128 =head2 accessors ref_allele()
 129
 130 reference allele nucleotide
 131
 132 =cut
 133
 134 has 'ref_allele' => ( isa => 'Str',
 135                    is  => 'rw',
 136     );
 137
 138 =head2 accessors alt_allele()
 139
 140 alternative allele nucleotide
 141
 142 =cut
 143
 144 has 'alt_allele' => ( isa => 'Str',
 145                    is  => 'rw',
 146     );
 147
 148
 149
 150 has 'qual' => ( isa => 'Str',
 151                 is  => 'rw',
 152     );
 153
 154 has 'filter' => ( isa => 'Str',
 155                   is => 'rw',
 156     );
 157
 158 has 'info' => ( isa => 'Str',
 159                 is => 'rw',
 160     );
 161
 162 has 'format' => ( isa => 'Str',
 163                   is => 'rw',
 164     );
 165
 166 has 'pAA' => ( isa => 'Num',
 167                is => 'rw',
 168               );
 169
 170 has 'pAB' => ( isa => 'Num',
 171                is => 'rw',
 172     );
 173
 174 has 'pBB' => ( isa => 'Num',
 175                is => 'rw',
 176     );
 177
 178 has 'raw' => ( isa => 'Str',
 179                is => 'rw',
 180     );
 181
 182 has 'error_probability' => ( isa => 'Num',
 183                              is => 'rw',
 184                              default => 0.001,
 185     );
 186
 187 =head2 get_score
 188
 189  Usage:        $ms->get_score('XYZ');
 190  Desc:         gets the marker score associated with XYZ
 191  Ret:
 192  Args:
 193  Side Effects:
 194  Example:
 195
 196 =cut
 197
 198 sub get_score {
 199     my $self = shift;
 200     my $accession = shift;
 201
 202     if (exists($self->scores()->{$accession})) {
 203         return $self->scores->{$accession};
 204     }
 205     else {
 206         warn "accession $accession has no associated score.\n";
 207         return undef;
 208     }
 209 }
 210
 211 sub snp_stats {
 212     my $self = shift;
 213
 214     my $good_snps = 0;
 215     my $invalid_snps = 0;
 216
 217     #print STDERR Dumper($self->snps());
 218     foreach my $a (@{$self->valid_accessions}) {
 219         my $snp = $self->snps()->{$a};
 220         if (!$snp) { print  "Warning! No snp found for accession $a!!!!\n"; }
 221         if ($snp && $snp->good_call()) {
 222             $good_snps++;
 223         }
 224         else {
 225             $invalid_snps++;
 226         }
 227     }
 228     return ($good_snps, $invalid_snps);
 229 }
 230
 231 sub calculate_allele_frequency_using_counts {
 232     my $self = shift;
 233
 234     my $total_c1 = 0;
 235     my $total_c2 = 0;
 236
 237     my $missing = 0;
 238     foreach my $k (@{$self->valid_accessions()}) {
 239         my $s = $self->snps->{$k};
 240         if (!$s) {
 241             print STDERR "Accession $k has no valid SNPs associated with it. Skipping.\n";
 242             $missing++;
 243             next();
 244         }
 245         $total_c1 += $s->ref_count();
 246         $total_c2 += $s->alt_count();
 247     }
 248
 249     if ($total_c1 + $total_c2 == 0) {
 250         return undef;
 251     }
 252
 253     my $allele_freq = $total_c1 / ($total_c1 + $total_c2);
 254
 255     my $pAA = $allele_freq **2;
 256     my $pAB = $allele_freq * (1 - $allele_freq) * 2 ;
 257     my $pBB = (1 - $allele_freq) **2;
 258
 259     $self->allele_freq($allele_freq);
 260     $self->pAA($pAA);
 261     $self->pAB($pAB);
 262     $self->pBB($pBB);
 263     #print STDERR "$missing accessions had no valid snps associated with them for snp ".$self->id()."\n";
 264     return $allele_freq;
 265 }
 266
 267  sub calculate_dosages {
 268      my $self = shift;
 269
 270      foreach my $k (keys %{$self->snps()}) {
 271         my $s = $self->snps()->{$k};
 272         $self->calculate_snp_dosage($s);
 273         my $dosage = $s->dosage() || "no dosage";
 274         #print STDERR "Dosage for SNP ".$s->id()."=$dosage for accession ".$s->accession()."\n";
 275      }
 276      #print STDERR Dumper($self->snps());
 277  }
 278
 279 =head2 function calculate_snp_dosage()
 280
 281 $snps->calculate_snp_dosage($snp)
 282 Calculates the SNP dosage of SNP $snp
 283 $snp is a CXGN::Genotype::SNP object
 284
 285 =cut
 286
 287 sub calculate_snp_dosage {
 288     my $self = shift;
 289     my $snp = shift;
 290     my $strict_dosage_filter = shift || 0;
 291
 292     my $error_probability = $self->error_probability();
 293
 294     #print STDERR "Working with strict dosage filter ($strict_dosage_filter) and error probability $error_probability\n";
 295     my $c1 = $snp->ref_count();
 296     my $c2 = $snp->alt_count();
 297
 298     #print STDERR "counts: $c1, $c2\n";
 299
 300     my $n = $c1 + $c2;
 301
 302     my $N1 = Math::BigInt->new($n);
 303     my $N2 = Math::BigInt->new($n);
 304
 305  #   print STDERR "$N1 bnok $c1 is: ". $N1->bnok($c1)."\n";
 306
 307     my $Nbnokc1 = $N1->bnok($c1)->numify();
 308     my $Nbnokc2 = $N2->bnok($c2)->numify();
 309
 310 #    print STDERR "NBnokc1: $Nbnokc1, NBnokc2 $Nbnokc2\n";
 311
 312     my $pDAA = $Nbnokc1 * ((1-$error_probability) ** $c1) * ($error_probability ** $c2);
 313     my $pDAB = $Nbnokc1 * (0.5 ** $c1) * (0.5 ** $c2);
 314     my $pDBB = $Nbnokc2 * ((1-$error_probability) ** $c2) * ($error_probability ** $c1);
 315
 316 #   print STDERR "pDAA: $pDAA pDAB $pDAB, pDBB $pDBB\n";
 317
 318     my $pSAA = $pDAA * ($self->pAA || 0);
 319     my $pSAB = $pDAB * ($self->pAB || 0);
 320     my $pSBB = $pDBB * ($self->pBB || 0);
 321
 322     if ($pSAA + $pSAB + $pSBB == 0) {
 323         return "NA";
 324     }
 325
 326     my $x = 1 / ($pSAA + $pSAB + $pSBB);
 327
 328     my $dosage = ($pSAB  + 2 * $pSBB) * $x;
 329
 330     if ($strict_dosage_filter) {
 331         if ( ($dosage > 0.1) && ($dosage < 0.9) ) {
 332             $dosage = "NA";
 333         }
 334         if ( ($dosage > 1.1) && ($dosage < 1.9) ) {
 335             $dosage = "NA";
 336         }
 337     }
 338
 339     $snp->dosage($dosage);
 340
 341     return $dosage;
 342 }
 343
 344 =head2 function hardy_weinberg_filter
 345
 346  $snps->hardy_weinberg_filter();
 347  returns a hash with the following keys:
 348    monomorphic   - 1 if the SNP is monomorphic
 349    allele_freq   - the allele frequency derived from counts
 350    chi           - the chi square value for hardy weinberg distribution
 351    scored_marker_fraction - the fraction of markers that were successfully scored
 352    heterozygote_count - the number of heterozygote SNPs in the set
 353
 354 =cut
 355
 356
 357 sub hardy_weinberg_filter {
 358     my $self = shift;
 359     my $dosages = shift; # ignored clones already removed
 360
 361     my %classes = ( AA => 0, AB => 0, BB => 0, NA => 0);
 362
 363     foreach my $d (@$dosages) {
 364         if (! defined($d) || $d eq "NA") {
 365             $classes{NA}++;
 366         }
 367         elsif ( ($d >= 0) && ($d <= 0.1) ) {
 368             $classes{AA}++;
 369         }
 370         elsif ( ($d >=0.9) && ($d <= 1.1) ) {
 371             $classes{AB}++;
 372         }
 373
 374         elsif (($d >=1.9) && ($d <= 2.0)) {
 375             $classes{BB}++;
 376         }
 377         else {
 378             #print STDERR "Dosage outlier: $d\n";
 379         }
 380
 381     }
 382
 383     print STDERR "Class counts: AA: $classes{AA}, BB: $classes{BB}, AB: $classes{AB}, NA: $classes{NA}\n";
 384
 385     if ( ( ($classes{AA} ==0) && ($classes{AB} ==0)) ||
 386         ( ($classes{BB} == 0) && ($classes{AB} ==0)) ) {
 387         return ( monomorphic => 1);
 388     }
 389
 390     my $total = $classes{AA} + $classes{AB} + $classes{BB};
 391
 392     my %score = ();
 393
 394     $score{total} = $total;
 395
 396     $score{scored_marker_fraction} = $total / (@$dosages);
 397
 398     #print STDERR "AA  $classes{AA}, AB $classes{AB}, BB $classes{BB} Total: $total\n";
 399     my $allele_freq = (2 * $classes{AA} + $classes{AB}) / (2 * $total);
 400
 401     $score{heterozygote_count} = $classes{AB};
 402
 403     $score{allele_freq} = $allele_freq;
 404
 405     my $expected = $allele_freq * (1-$allele_freq) * 2 * $total;
 406
 407     #print STDERR "TOTAL: $total\n";
 408     my $x = ($classes{AB} - $expected)**2 / $expected;
 409
 410     # only do the chi square test if the number of heterozygotes is larger than expected
 411     #
 412     if ($classes{AB} > $expected) {
 413         $score{chi} = $x;
 414     }
 415     else {
 416         $score{chi} = 0;
 417     }
 418
 419     return %score;
 420 }
 421
 422
 423 __PACKAGE__->meta->make_immutable;
 424
 425
 426 1;